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Abstract: We determine the MS charm quark mass from a charmonium QCD sum 
rules analysis. On the theoretical side we use input from perturbation theory at Q(a%). 
Improvements with respect to previous 0{o^) analyses include (1) an account of all 



available e + e~ hadronic cross section data and (2) a thorough analysis of perturba- 
tive uncertainties. Using a data clustering method to combine hadronic cross section 
data sets from different measurements we demonstrate that using all available exper- 
imental data up to cm. energies of 10.538 GeV allows for determinations of experi- 
mental moments and their correlations with small errors and that there is no need to 
rely on theoretical input above the charmonium resonances. We also show that good 
convergence properties of the perturbative series for the theoretical sum rule moments 
need to be considered with some care when extracting the charm mass and demonstrate 
how to set up an appropriate set of scheme variations to obtain an appropriate esti- 
mate of the perturbative uncertainty. As the final outcome of our analysis we obtain 
m c (m c ) = 1.277 ± 0.006 sta t ± 0.013 syst ± 0.019p Crt ± 0.009 Qs ± 0.002 (GG) GeV. The pertur- 
bative error is an order of magnitude larger than the one obtained in previous O(a^) sum 
rule analyses. 
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1. Introduction 

Accurate determinations of the charm quark mass are an important ingredient in the 
prediction of inclusive and radiative B decays or exclusive kaon decays such as K — > irvV. 
Since these decays are instruments to either measure CKM matrix elements or to search 
for new physics effects, appropriate and realistic estimates of the uncertainties are also an 
important element of these analyses Q. 

One of the most powerful methods to determine the charm quark mass is based on sum 
rules for the charm-anticharm production rate in e + e~ annihilation ||. Here, moments of 
the correlation function of two charm vector currents at zero momentum transfer 

{9^q 2 -q„qu) I% 2 ) = -i J dx e iqx ( \T j^x)j v {Q)\ ) , 
f(x) = i){x)^i){x) , 
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Q c being the charm quark electric charge, can be related to weighted integrals of the 
normalized charm cross section 



Mn = I 4^T R e+e-^cc+x(s), (1.2) 



Re+e-^cc+x{ s ) 



• cc+X 



00 



which can be obtained from experiments. For small values of n such that m c /n > Aqcd 
the theoretical moments M* h can be computed in an operator product expansion (OPE) 
where the dominant part is provided by perturbative QCD supplemented by small vacuum 
condensates that parametrize nonperturbative effects (3|, |4j. The leading gluon condensate 
power correction term has a surprisingly small numerical effect and is essentially negligible 
for the numerical analysis as long as n is small. 

This allows to determine the charm mass in a short distance scheme such as MS 
to high precision. This method to determine the MS charm mass is frequently called 
charmonium sum rules. For the theoretical moments the perturbative part of the OPE 
is known at 0(a®) and 0(a s ) for any value of n ||. At 0{a 2 s ) the first 30 moments 
are known §, 0, and to 0(a 3 s ) for n = 1 @, @, n = 2 Q, and n = 3 0. Higher 
moments at O(a^) have been determined by a semianalytical procedure [12, O] (see also 



|l4|]). The Wilson coefficient of the gluon condensate contribution is known to 0(a s ) 
|15|| . On the experimental side the total hadronic cross section in e + e~ annihilation is 
known from various experimental measurements for cm. energies up to 10.538 GeV. None 
of the experimental analyses actually ranges over the entire energy region between the 
charmonium region and 10.538 GeV, but different analyses overlapping in energy exist 



such that energies up to 10.538 GeV are completely covered |16|, [17|, U, f[9[ |2Cj, |21, ^2 



H ||, ||, ||, |7|, H, ||, |(| |l], ||, H. 1 Interestingly, to the best of our knowledge, the 
complete set of all available experimental data on the hadronic cross section has never been 
used in previous charmonium sum rule analyses to determine the experimental moments. 
Rather, sum rule analyses have relied heavily on theoretical input using different approaches 
to determine the corresponding "experimental error" and intrinsically leading to a sizable 
modeling uncertainty for energy regions below 10.538 GeV for low values of n [p 7 ^]. 

The most recent charmonium sum rule analysis based on Eqs. ( |l . l| ) and (|l.2| ), carried 



out by Kiihn et al. [35, f3(| using input from perturbative QCD (pQCD) at O(a^) for the 
perturbative contribution, obtained m c (m c ) = 1279± (2) per t± (9) ex p± (9) Qs ± (1)(gg> MeV 
where the first error is the perturbative uncertainty and the second is the experimental one. 
The third and the fourth uncertainties come from a s and the gluon condensate correction, 
respectively. To our knowledge this result, the outcome of similar analyses in Ref. B and 



*As a word of caution we mention that, except for the contributions of the J/ip and ip' resonances, no 
experimental separation of the charm and non-charm contributions in the hadronic cross section has been 
provided in available data, although charm-tagged measurements are possible, see e.g. Ref. (CLEO 
collaboration). So the charm pair production rate from above the J /if) and tp' that enters the charmonium 
sum rules in Eq. jl.2| ) is usually obtained partly from the measured total i?-ratio with theory motivated 
subtractions of the non-charm rate, and partly by using theory predictions for the charm production rate. 
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by Boughezal, Czakon and Schutzmeier || 2 , and a closely related analysis based on lattice 
results instead of data for pseudoscalar moments ^] represent the analyses with the 
highest precision achieved so far in the literature. If confirmed, any further investigations 
and attempts concerning a more precise charm quark MS mass would likely be irrelevant 
for any foreseeable future. 

We therefore find it warranted to reexamine the charmonium sum rule analysis with 
special attention on the way how perturbative and experimental uncertainties have been 
treated in Refs. |?5], [36|. A closer look into their analysis reveals that the quoted pertur- 
bative uncertainty results from a specific way to arrange the a s expansion for the charm 
mass extractions and, in addition, by setting the MS renormalization scales in a s and in 
the charm mass (which we call /x Q and fi m , respectively) equal to each other (i.e., they use 
[J>a = fJ"m)- Moreover, concerning the experimental moments, only data up to y/s = 4.8 GeV 
from the BES experiments 17, |2(| were used, while for y/s > 4.8 GeV perturbative QCD 



predictions were employed. Conceptually this approach is somewhat related to the method 
of finite energy sum rules (see e.g. Ref. pOp ), which we, however, do not discuss in this 
work. While this approach might be justified to estimate the overall nominal contribution 
for the experimental moments from y/s > 4.8 GeV, since perturbative QCD predictions 
describe quite well the measured total hadronic cross section outside the resonance regions, 
it is certainly not a suitable method to determine the experimental uncertainty. Since the 
region y/s > 4.8 GeV constitutes about 30% of the first moment Mi, which is theoretically 
most reliable, this approach contains a significant intrinsic model dependence that cannot 
be quantified unambiguously. 

In this work we reexamine the charmonium sum rules analysis for low values of n using 
the latest 0(a^) perturbative results, and we implement improvements which concern the 
two issues just mentioned: 

1. We analyze several different types of perturbative expansions and examine in detail 
how the result for the MS charm mass depends on independent choices of /i Q and 
fj, m . We show in particular that the interplay of the perturbative expansion and the 
scale setting /j, a = \i m used in previous 0{a\) analyses leads to sizable cancellations 
of the dependence on fi a and \i m that in the light of our new analysis has to be 
considered as accidental. As the outcome of our analysis we estimate the current 
0(a>g) perturbative error as around 20 MeV, which is an order of magnitude larger 



than that of Refs. |35|, 36 



Using a clustering method [41, 42, 43] to combine correlated data from many differ- 
ent experimental measurements we show that the e + e~ total hadronic cross section 
relevant for the charmonium sum rules can be determined with a complete coverage 
of center of mass energies above the J/ip and tp' resonances up to 10.538 GeV. Con- 
servatively estimated modeling uncertainties coming from the energy range above 



2 Since the analyses of Refs. [g, |9| were based on outdated and less precise data for the J/ip and ip' 
electronic partial widths |$7||, we frequently only compare our numerical results with those of Refs. pr| . 
However the perturbative input of Refs. |35[ |3^] and ours is identical. 
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10.538 GeV then only lead to an insignificant contribution to the total uncertainty of 
the experimental moments. 

This paper is organized as follows: In Sec. ^ we introduce the theoretical framework and 
review the current status of perturbative computations. We also show various equivalent 
ways of arranging the perturbative series in a s for the charm mass. Finally we discuss how 
to properly estimate theoretical uncertainties due to the truncation of the perturbative 
series. In Sec. |3] we present all the experimental information that goes into our analysis. 
We discuss a clustering fit procedure that allows to combine data from different experiments 
accounting for their correlation and show the results. In Sec. Qwe carry out the numerical 
charm mass analysis concentrating on the first moment using arbitrary values of a s , and 
we present our final charm mass result. In Appendix A we present more details on the 
outcome of our clustering fit procedure for the charm i?-ratio, and in Appendix B we 
prove the equivalence of different versions of x 2 functions when auxiliary fit parameters are 
employed. 



2. Theoretical Input 

2.1 Perturbative Contribution 

The moments of the vector current correlator are defined in Eq. ( |l.ip . Their perturbative 
contribution in the framework of the OPE has a non-linear dependence on the charm 
quark mass. Thus in principle no conceptual preference can be imposed on any of the 
possible perturbative series that arises when solving for the charm mass. As a consequence, 
different versions of the expansion should be considered to obtain reliable estimates of 
the perturbative uncertainty. As indicated in Sec. |l| we use in the following fj, a as the 
renormalization scale in a s and fi m as the renormalization scale in the MS charm quark 
mass m c . 



(a) Standard fixed-order expansion 

Writing the perturbative vacuum polarization function as 

1 oo 

II pert (9 2 } a a (Ma),?ScW^a,/im) = J^2Q2 E > ( 2A ) 

c n=0 

we have for the perturbative moments M^ ert 

pert = i v ( 1 a* in* ( m 6 ( . (2.2) 

This is the standard fixed-order expression for the perturbative moments. At O(a^) the 
coefficients C^'g were recently determined forn = lpl0],n = 2 |l(| , n = 3 JT1J and higher 
12, |l3|, U]. We refer to Ref. [||, 0] for the coefficients at 0{a 2 s ). We have summarized the 



numerical expressions for the C°'j coefficients of the first four moments in Tab. |l]. 



The standard fixed-order expansion in Eq. (2^2) is the common way to represent the 



perturbative moments. However, written in this form the non-linear dependence on m c 
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does for some values of the experimental moments and the renormalization scales not yield 
numerical solutions for m r . 3 
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Table 1: Numerical values of the coefficients for Eq. (2.2). (Standard fixed-order expansion). 



(b) Linearized expansion 

Concerning the charm mass dependence, a more linear way to organize the perturbative 
expansion is to take the 2n-th root of Eq. (2.2): 



M^V /2n = ) V ( ^Yc^ln" ( ha* ( , (2.3) 



i.a.b 



K J V Mm / V M 



or equivalently 

2 ^ M th,pert\ . ab V 



The coefficients C®\ using again jx a for the renormalization scale in a s and fi m for the 



renormalization scale in the MS charm mass are given in Tab. g. Although relation ( |2~j| ) 
involves a non-linear dependence on fn c , we find that it always has a numerical solution. 

(c) Iterative linearized expansion 



For the standard and the linearized expansions in Eqs. (2.2) and Q2.3]) one searches for 



3 This can to happen at any order for n > 1 and for (i m ~ 3 GeV. 
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Table 2: Numerical values of the coefficients for Eq. ( |2.3[ ). (Linearized expansion). 



numerical solutions of the charm mass m c ((i m ) keeping the exact mass dependence on the 
respective RHS at each order in a s . An alternative way that is consistent within pertur- 
bation theory is to solve for m c (^ m ) iteratively order- by-order supplementing appropriate 
lower order values for rn c (/j, m ) in higher order perturbative coefficients. To be more explicit, 
we describe the method in the following. As the basis for the iterative expansion carried 
out in our analysis we use the linearized expansion of Eq. fl2.4|) . 

In the first step we determine m c (// m ) employing the tree-level relation 



1 



M* h ' pert ) 1/2? 



C 



o.o 

n,0 ' 



(2.5) 



giving the tree-level charm mass mc°^(/i m ). In the next step one employs the relation 



2[Mn 



_J j ^0,0 , ttsiVa) 



\ H'm 



, (2.6) 



to determine the 0(a s ) charm mass m' 1 ' (/i m ). In the 0(a s ) terms on the RHS of Eq. 
the tree- level charm mass is used, which is consistent to 0(a s ). At 0(a 2 s ) for the 



determination of rrii\^ m ) one uses m^' for the 0(aj) coefficient and m^'f/im) for the 



,(0) 
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0(a s ) correction, which in the strict a s expansion yields 



ml 2 ) 
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In" 



(0)2 



Or, 



/'o 



(2.7) 



Here the second line contains the derivative of the 0(a s ) terms with respect to the charm 
mass. The determination of the 0{oP s ) charm mass rrffi (nm) is then carried out in an 
analogous way involving the second (first) derivative with respect to the mass in the 0(a s ) 
{0{a 2 s )) correction and using again mc°^(// m ) for the 0{o? s ) coefficient. 
In general we can write the iterative expansion as follows: 



i,a,b 



7T 



d a 6 ln a 



m 



(0)2^ 



ln b 



m. 



<0)2' 



(2.8) 



where the numerical value of the coefficients are collected in Tab. [3|. 

The iterative way to treat the perturbative series for the charm mass has the advantage 
that solving for the charm mass involves equations that are strictly linear in the charm 
mass at any order of the a s expansion and thus always have solutions. In this way any 
possible influence on the analysis arising from a non-linear dependence is eliminated. 



(d) Contour improved expansion 

For the expansion methods (a)-(c) the moments and the charm quark mass are computed 
for a fixed choice of the renormalization scale fi a in the strong coupling a s . In analogy to 
the contour improved methods used for r-decays (see e.g. Refs. || H, EH Hi) one 
can employ a path-dependent ix a in the contour integration that defines the perturbative 
moments |54j, see Fig. [I], 



^c,pcrt _ 67rQ e 



d.S 



T I% , m c )),m c (n m ),fi c a (s, m c ), /x. 



(2.9) 



Due to the independence of the moments on [i a and since no large logarithms are being 
generated anywhere for a path with distance of order fn c from the cut on the real axis, 
this method is a viable alternative to carry out the perturbative expansion. The different 
orders in the expansion of the contour improved moments M^' pert are generated from the 
fixed-order a s expansion of the vacuum polarization function II in Eq. (|2.9|) . 
A useful path-dependent choice for [i c a is given by [53 
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Table 3: Numerical values of the coefficients for Eq. (Ej). (Iterative linearized expansion). 
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Figure 1: Path of integration in the complex s-plane for the computation of the moments. 



which implements a modified weighting of threshold versus high energy contributions. It 
is straightforward to prove that the resulting moments Mn' pert can be obtained from the 
small-g 2 expansion of the perturbative vacuum polarization function using as the renor- 
malization scale of a s , 

oo 

^^(g^a s (^(g^m2)),m c ( Mm ),^( 9 2 ,m c 2 )^ m ) = £ q 2n M^ ert . (2.11) 

n=0 

From Eq. (2.11) we can see that the M^' port can be derived from the expressions for the 
fixed-order moments Mm rt with m < n given in Eq. (|2 . 2[) . 4 They also depend on the QCD 



4 This works in general as long as the path dependent /i^g 2 , m 2 .) does not produce a spurious cut in a s 
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/3-function and its derivatives, which arise in the small-g 2 expansion of a s (^(g 2 , m^)). 
Note that the /3-function that has to be employed must be exactly the same that is used 
for the contour integration in Eq. (|2.S|). Expanding the dependence of the /3-function on 
a s strictly in fixed-order one recovers the fixed-order moments M^ ert . So the dependence 
of the contour improved moments pert on the fixed-order moments Mm ert with m < n 
is only residual due to the truncation of the a s series, representing yet another alternative 
parametrization of higher order perturbative corrections. The contour improved moments 
do in particular have a residual dependence on the UV-subtraction scheme for the vacuum 
polarization function, i.e. on 11(0) = MQ Grt . Using the "on-shell" scheme with n(0) = one 



finds that Mi' peit 



Mf ert . For our analysis we employ the MS scheme for n(0) defined 



for /x = m c {m c ). Expressed in terms of a s (fi a ) and m c (/x m ) it has the form 
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The numerical values for the coefficients Cq'- can be found in Tab. ^. Using contour 
improved moments to determine m c also involves non-linear relations, which implies that 
in some cases there is no solution. Again this can happen for n => 1 and for fi m ~ 3 GeV 
at any order. 
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Table 4: Numerical values of the coefficients for Eq. (2.12). (11(0) in the MS scheme.) 



2.2 Gluon Condensate Contribution 



The dominant subleading contribution in the OPE for the theory moments M* h is arising 
from the gluon condensate. Writing 
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(2.13) 



and neglecting any higher order power-suppressed condensate contributions of the OPE we 
have |, H, HI 
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• (2-14) 



We have written Eq. ( p,14| ) using the renormalization group invariant (RGI) scheme for 
the gluon condensate InJ and the MS charm quark mass. In Ref. it is suggested 



starting at q 2 — and running towards — oo. This condition is implemented into Eq. (2.10) by (Ha) 2 being 
negative along the physical cut of the vacuum polarization function above the charm pair threshold. 
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that the Wilson coefficient of the gluon condensate should be expressed in terms of the 
pole rather than the MS mass based on the observation that the pole mass leads to a 
condensate correction that is numerically quite stable for higher moments. We disagree 
with this prescription due to the infrared renormalon issue in the pole mass and since it 
is expected that the condensate corrections for the n-th moment increases with n. On the 
other hand, for a sum rule analysis of the first moment with n = 1 this issue is absolutely 
irrelevant due to the small size of the gluon condensate contribution compared to the level 
of the remaining perturbative uncertainties. Numerical values for the coefficients a % n are 
given in Tab. || for n = 1, 2, 3, 4. For the RGI gluon condensate we adopt 53] 



—G 2 \ 

TT I RGI 



0.006 ± 0.012 GeV 4 



(2.15) 



The overall contribution of the gluon condensate correction in Eq. ( |2.14| ) in the charm 
quark mass analysis is quite small. Its contribution to the moments amounts to around 
0.2%, 0.6%, 2%, and 3% for the first four moments, respectively. For n = 1 it leads to a 
correction in the MS charm quark mass at the level of 1 MeV and is an order of magnitude 
smaller than our perturbative uncertainty. We therefore ignore the condensates correction 



for the discussion of the perturbative uncertainties in Sec. |2.4j . Its contribution is, however, 
included in the final charm mass results presented in Sec. 0. 



0,0 



1,0 

a n 

1,1 
a n 



n 



1 



n 



n 



n 



-16.042 

-15.0278 

-96.2521 



-26.7367 
13.0055 
-213.894 



-38.8898 
78.710 
-388.898 



-52.3516 
190.935 
-628.219 



Table 5: Numerical values for the coefficients of Eq. (2.14). (Gluon condensate contribution). 



2.3 Running Coupling and Mass 

The analysis of the charmonium sum rules naturally involves renormalization scales around 
the charm mass, \i ~ m c ~ 1.3 GeV, which are close to the limits of a perturbative treat- 
ment. In fact, parametrically, the typical scale relevant for the perturbative computation 
of the n-th moment M n h is of order u ~ m c /n (see e.g. Ref. [54]) because the energy range 
of the smearing associated to the weight function l/s n+1 in Eq. (1.2) decreases with n. We 
will therefore use n = 1 for our final numerical analysis. Moreover, it is common practice 
to quote the MS charm mass m c (m c ), i.e. for the scale choice [i m = m c . It is therefore 
useful to have a look at the quality of the perturbative behavior of the renormalization 
group evolution of the strong MS coupling a s and the MS charm quark mass. 

In Fig. [2(a)p we have displayed a^ 3LL (^)/a^ kLL (u) using for a^ kLL the (k + l)-loop 



QCD /3-function and the respective exact numerical solution for a s (3 GeV) = 0.2535 as the 
common reference point. We see that the convergence of the lower order results towards the 



5 In this examination and throughout our other analyses we use the MS renormalization group equations 
with nf = 4 active running flavors. 
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4-loop evolution is very good even down to scales of around 1 GeV. The curves indicate that 
the remaining relative perturbative uncertainty in the 4-loop evolution might be substan- 
tially smaller than 1% for scales down to rn c ~ 1.3 GeV. It is also instructive to examine 
the evolution using a fixed-order expansion. In Fig. gEJ we display af LL (n)/a i J n) (fJ-) 



where ai"^(/x) is the C(a™ +1 ) fixed-order expression for a s (/i) using the reference value 
a s (3 GeV) = 0.2535 as the expansion parameter. The convergence of the fixed-order ex- 
pansion for a s (fi) towards the exact N 3 LL numerical solution a^ 3LL (/i) is somewhat worse 
compared to the renormalization group resummed results since the deviation of the ratio 
from one is in general larger. However, convergence is clearly visible. In particular there 
are not any signs of instabilities. It therefore seems to be safe to use renormalization scales 
down to the charm mass and associated renormalization scale variations as an instrument 
to estimate the perturbative uncertainties. 



In Figs. 3(a) and 3(b) an analogous analysis has been carried out for the MS charm 



quark mass. In Fig. 3(a)| m^ iLL (^)/m^ LL (//) is plotted for k = 0, 1, 2, 3 using the exact nu- 



merical solutions of the (k + l)-loop renormalization group equations and rn c (/x = 3 GeV) 
as the respective reference value. 6 Compared to the Fig. ^ we observe a very similar con- 
vergence. In Fig. 3(b)| , finally, we show LL (|ii)M[ m '(/i), where mi m ^(/x) is the 0{a™) 



fixed-order expression for m c (//) using a s (3 GeV) = 0.2535 as the expansion parameter. 
Again, the convergence towards the exact N 3 LL evolved result is very similar to the cor- 
responding results for the strong coupling, and we find again no evidence for perturbative 
instabilities. Of course the corrections are somewhat larger when the fixed-order expan- 
sion is employed. We therefore conclude that perturbative evolution and renormalization 
scale variations for the MS charm quark mass can be safely used down to scales above 
m c ~ 1.3 GeV. One should of course mention that the lines in Fig. 3(b) also give an indi- 



cation about the expected size of scale variations depending on the range of the variations. 
Scales above 2 GeV can lead to sub-MeV variations, while scales down to the charm mass 
will result at best in percent precision (i.e. O(10 MeV)). 

2.4 Perturbative Uncertainties in the MS Charm Mass 

In this section we discuss in detail the perturbative series for the determination of the MS 
charm mass m c and how to set up an adequate scale variation to estimate the perturbative 
uncertainty. In the previous subsections we have discussed four different ways to carry out 
the perturbative expansion and we presented the corresponding order-by-order analytic 
expressions. As described there, we can determine at each order of the perturbative expan- 
sion for the moments M^ 1 a value for rn c (/i m ) which also has a residual dependence on 
the renormalization scale used for a s . To compare the different mass determinations we 
then evolve m c (/i m ) to obtain m c (m c ) using the 4-loop renormalization group equations for 



the mass and the strong coupling |55|, [56L 57, 58, p9|. The obtained value of m c {m c ) thus 



6 The numerical value of m c (/i = 3 GeV) is actually irrelevant, since the running involves m c only in a 
linear way and exactly cancels in the ratio. 

7 For the discussions in this section we use ai" f 5 \mz) = 0.118 (ai nf 4 '(4.2 GeV) = 0.2245) as an 
input using five-to-four flavor matching at 4.2 GeV. For the first moment employed for the charm mass fits 
we use Mi = 0.2138 GeV -1 , corresponding to the standard selection of data sets, as described in Sec. |^. 
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1.00 
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3.0 



(a) 



(b) 



Figure 2: Results for a™ 3jjh (//) / a^ hh (jjl) (a) and af hL (n)/a ( s k) (/i) (b), where o^ kLL stands for 
the (fc + l)-loop running coupling constant and ai*^ is the corresponding 0{pts +1 ) fixed-order 
expression for a s . All orders are run from the common point a s (3 GeV) = 0.2535. 
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Figure 3: Results for mf LL (,u)/mf LL (/i) (a) and mf* LL (ii)/m^ fe) (/x) (b), where to^ Kll stands 



N 3 LL / 



1.5 



2.0 2.5 



(b) 



3.0 



for the [k + l)-loop running MS charm mass and mi is the 0(a 



fixed-order expression. 



has a residual dependence on the scales fj, m and fj, a , on the order of perturbation theory 
and on the expansion method. 8 For the results we can therefore use the notation 



m c (m c )[ 

where n = 0, 1, 2, 3 indicates perturbation theory at 0(a™) and 

a (fixed-order expansion), 

b (linearized expansion), 

c (iterative expansion), 

d (contour improved expansion). 



(2.16) 



(2.17) 



To initiate the discussion, we show in Fig. ^results for m c (m c ) at 0(a™) for expansions 
a-d using fj, m = fj, a (upper four graphs) and using \x m = m c (m c ) (lower four graphs). For 



8 Of course the extracted mass depends on the moment considered as well. Since in our analysis we focus 
on the first moment only, for simplicity we drop that label. 
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m c {m c ) 



1.35 



1.30 



1.25 



1.20 




0(<) 

Figure 4: Results for m c (m c ) at various orders, for methods a (graphs 1 and 5), b (2,6), c (3,7), 
and d (4,8), setting \x a = (j, m (graphs 1-4) and setting /x m = m c (m c ) (5-8). The shaded regions 
arise from the variation 2 GeV < fi a < 4 GcV. 



1.30p 

1.29 

1.28 



m c (m c ) various methods a(mz)=0.118 



I 



m c {m c ) 



2 GeV</j„<4 GeV 

/J,„=3 GeV 

2 GeV<// ( ,<4 GeV 

3 

fi m =m c {m c ) 
2GeV<^ a <4 GeV 

m c (m c )<fi m <4 GeV 
m c (m c )<fi a <4 GeV 



1.27 

1.26 Fixed Order 
1.25^ 



II 



Contour 
Improved 



Linearized 



Iterative 



Method 

Figure 5: Estimates of the perturbativc error at 0(a^). We show the correlated [i m = /i Q 
variation, (orange), setting fi m — m c (m c ) and setting /x m = 3 GeV, (blue and red, respectively), 
and the double scale variation (magenta). 



each method and order we have displayed the range of rn c (jn c ) values for a variation of 
2 GeV < fi a < 4 GeV, which corresponds to the scale variation employed in Refs. [||. |{J 
35, 36 [. Their analysis used the fixed-order expansion with the setting fjj m — and is 
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represented by graph 1. We make several observations: 



(i) Choosing \x m and fi a both larger than 2 GeV makes the m c (m c ) value decrease with 
the order of perturbation theory. 

(ii) Choosing fi m smaller than 1.5 GeV and \i a larger than 2 GeV makes the m c (m c ) 
value increase with the order of perturbation theory. 

(hi) For most choices of the scale setting and the expansion method the spread of the 
^c(^c) values from the variation of fi a does not decrease in any substantial way 
with the order. However, viewing all methods and scale setting choices collectively a 
very good convergence is observed. 

We have checked that these statements apply also in general beyond the specific cases 
displayed in Fig. ||. 

Quite conspicuous results are obtained for the scale choice [i m = fi a for the fixed- 
order (graph 1) and linearized expansions (graph 2). Here, extremely small variations in 
m c (m c ) are obtained. They amount to 1.8 MeV (4 MeV) and 0.6 MeV (1.4 MeV) at 
order a 2 s and a^, respectively, for the fixed-order expansions (linearized expansions). We 
note that our scale variation for the fixed-order expansion at O(a^) is consistent with 
the corresponding numbers quoted in Ref. [||, 9, where the 0(a^) corrections to the first 
moment were computed, 9 but differs from the scale variations given in Refs. [35, 36] which 
also quoted numerical results from the fixed-order expansion. Interestingly the 0{a 2 s ) 
and 0(aP s ) variations we find do not overlap and the O(a^) ranges appear to be highly 
inconsistent with the O(cef) results from the iterative method (graph 3). A visual display 
of scale variations obtained from the four expansion methods at 0(a^) with different types 
of variation methods is given in Fig. ||. 

An illustrative way to demonstrate how a small scale variation can arise is given in 
Fig. 6(a) - 6(d). For all four expansion methods contour curves of constant m c (m c ) are 



displayed as a function of (the residual dependence on) fj, m and fj, a . For the fixed-order 
(a) and the linearized expansions (b) we see that there are contour lines closely along the 
diagonal fi m = [i a . For the fixed-order expansion (a) this feature is almost exact and thus 
explains the extremely small scale-dependence seen in graph 1 of Fig. |5|. For the linearized 
expansion (b) this feature is somewhat less exact and reflected in the slightly larger scale 
variation seen in graph 2 of Fig. ^. On the other hand, the contour lines of the iterative 
(c) and contour improved (d) expansions have a large angle with respect to the diagonal 
Mm — leading to the much larger scale variations of 9 and 5 MeV, respectively, at O(a^) 
visible in graphs 3 and 4 of Fig. ||. The contour plots shown in Fig. || also show the variation 
along the line fj, m or fj, a = m c (rn c ) (border of gray shaded areas). Here, the contour lines 
are relatively dense leading to scale variation of around 15 MeV at 0{a? s ). 
Overall, we draw the following conclusions: 



9 We are grateful to Thomas Schutzmeier for confirming agreement of the results of our O(a'l) fixed-order 
code with theirs of Ref. [HI . 
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Mm 
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Figure 6: Contour plots for m c (m c ) as a function of fi a and ^ m at 0(ctg), for methods (a)-(d). 
The shaded areas represent regions with /i m , \i a < m c (m c ), and are excluded of our analysis. 

(1) The small scale variations observed for the fixed-order and linearized expansions for 
A*m = /■*£* result from strong cancellations of the individual fj, m and (i a dependences 
that arise for this correlation. 

(2) Other correlations between fx m and fx a that do not generate large logarithms do not 
lead to such cancellations. One therefore has to consider the small scale variations 
observed for \i m = \i a in the fixed-order and linearized expansions as accidental. 



(3) For an adequate estimate of perturbative uncertainties specific correlations between 
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fj, m and (jL a that are along contour lines of constant m c (m c ) have to be avoided. 
Moreover, adequate independent variations of \i m and fj, a should not induce large 
logarithms. 

As the outcome of this discussion we adopt for our charm mass analysis an independent 
and uncorrelated variation of fj, m and fi a in the range 

rn c (m c ) < Hm, Ha < 4 GeV , ( 2 -18) 

in order to estimate perturbative uncertainties. The excluded region fj, m ,pi a < m c (fn c ) in 
the Hm-fJ-a plane in the contour plots of Fig. |(| is indicated by the gray shaded areas. This 
two-dimensional variation avoids accidental cancellations from correlated jU m -/U Q variations, 
large logarithms involving ratios of the scales fn c (fj, m ), fj, m , [i a and remains well in the 
validity ranges for the perturbative renormalization group evolution of the MS charm quark 
mass and a s (see Sec. |2,3| ). The range of Eq. fl2.18|) is also consistent with a scale variation 
[i ~ 2m c x (1/2,1,2) one might consider as the standard choice with respect to the to 
particle threshold located at \fs = 2fn c . As a comparison the range // = (3 ± 1) GeV 
corresponds to fj, = 2m c (0.8, 1.2, 1.6). We emphasize that the lower boundary m c (m c ) ~ 
1.25 GeV is also reasonable as it represents the common flavor matching scale where gauge 
coupling evolution remains smooth up to NLL order. We do not see any evidence for 
perturbation theory for the first moment M\ being unstable at the charm mass scale, in 



contrast to claims made in Ref. [36]. 



In Fig. [7(a) we show the ranges of m c (m c ) at 0(al' 2,3 ) for the four expansion methods 
employing the scale variations of Eq. fl2.18|) . We see that all four expansion methods now 
lead to equivalent results. The variations are also compatible with the overall variations 
shown in Fig. |j. At O(a^) we obtain a scale variation for m c (rn c ) of around 20 MeV. This 
is an order of magnitude larger than the perturbative uncertainties quoted in Refs. []|, || 



35, pq| . In Fig. 7(b) we have displayed the corresponding results for m c (3 GeV). At 0{a 



3\ 



they also exhibit a scale variation of around 20 MeV. 

3. Experimental Data 
3.1 Data Collections 
Narrow resonances 

Below the open charm threshold there are the J /ip and tp' narrow charmonium resonances. 



Their masses, and electronic widths are taken from the PDG [60] and are collected in Tab. y 
together with the value of the MS QED coupling at their masses. The total widths are 
not relevant since we use the narrow width approximation for their contributions to the 
moments. The uncertainty for the contribution to the moments coming from the masses 
can be neglected. 

Threshold and data continuum region 

The open charm threshold is located at yfs = 3.73 GeV. We call the energies from just 
below the threshold and up to 5 GeV the threshold region, and the region between 5 GeV 
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O(o*) 



(a) 
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Figure 7: Results for m c (m c ) (a) and m c (3 GeV) (b), for methods a-d, at orders C(a^ 2,3 ). At 
0{a\) the error bars for fixed-order method extend down to 0.93 GeV and 0.62 GeV, for m c (m c ) 
and m c (3 GeV), respectively. 







M (GeV) 
Tee (keV) 
(a/a{M)) 2 


3.096916(11) 3.686093(34) 
5.55(14) 2.48(6) 
0.957785 0.95554 



Table 6: Masses and electronic widths |30[ of the narrow charmonium resonances and effective 
electromagnetic coupling |$(|. a = 1/137.035999084(51) stands for the fine structure constant, and 
a(M) stands for the MS electromagnetic coupling at the scale M. 



and 10.538 GeV, where the production rate is dominated by multiparticle final states 
the data continuum region. In these regions quite a variety of measurements of the to- 
tal hadronic cross section exist from BES |6|, [17J, |l|, 0, fl], CrystalBall |22], g3|, 

" 13, Ba||. Taken 



CLEO gj, H, H, |7|, MDl ||], PLUTO |g|, and MARKI and II |fj, [M 



together, the entire energy region up to 10.538 GeV is densely covered with total cross sec- 
tion measurements from these 19 data sets. 10 The measurements from BES and CLEO have 
the smallest uncertainties. They do, however, not cover the region between 5 and 7 GeV. 
Here CrystalBall and MARKI and II have contributed measurements albeit with somewhat 
larger uncertainties. The statistical and total systematical uncertainties of the measure- 
ments can be extracted from the respective publications. For some data sets the amount of 
uncorrelated and correlated systematical uncertainties is given separately (BES (17], [l!| ^Cj . 
CrystalBall pi, 



], CLEO p§, MARKI and II [33, 32], MDl pa]) while for all the other 



data sets only combined systematical uncertainties are quoted. All these data sets are 
shown in Figs. |8[ where the displayed error bars represent the (quadratically) combined 
statistical and systematical uncertainties. 

Interestingly, none of the previous charm mass analyses, to the best of our knowledge, 



10 There are 18 references quoted since Ref. |2^] provides results from two independent runs that we treat 
as two different data sets. 
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Figure 8: Experimental data. In the two top figures we show BES and CLEO data in the non- 
charm, low charm and threshold regions (a) and low charm region and threshold region (b). The 
two pictures in the middle show less precise data in the threshold region: CrystalBall and MARKII 
(c) and MARKI and PLUTO (d). The two pictures on the bottom show data in the continuum 
region: below 8 GeV data from CrystalBall, MARKI and PLUTO (e) and above 8 GeV data from 
CLEO, MD1 and PLUTO (f). 



ever used the complete set of available data. As examples, Bodenstein et al. [61] used 



data sets [16, 17, 20, £7| from BES and CLEO. Jamin and Hoang [34] used the data sets 



of Refs. 0, HI and H from BES, MD1 and CLEO, covering the regions 2 GeV < E < 
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Numbering 


Reference 


Experiment 


Year 


Systematical 


Data points 


1 


in 


BES 


2000 


no splitting given 


6 


2 




17] 


BES 


2002 


splitting given 


85 


3 


in 


BES 


2004 


only one point 


1 


4 


@ 


BES 


2006 


splitting given 


3 


5 


@ 


BES 


2006 


splitting given 


68 


6 




21] 


BES 


2009 


no splitting given 


3 


7 


m 


CrystalBall 


1986 


splitting given 


98 


8 


m 


CrystalBall (Run 1) 


1990 


splitting given 


4 


9 


m 


CrystalBall 


1990 


splitting given 


11 


10 




24] 


CLEO 


1979 


only one point 


1 


11 


m 


CLEO 


1998 


only one point 


1 


12 


in 


CLEO 


2007 


splitting given 


7 


13 




27] 


CLEO 


2008 


no splitting given 


13 


14 


m 


MD1 


1996 


splitting given 


31 


15 


m 


PLUTO 


1982 


no splitting given 


45 


16 


13 


MARKI 


1976 


no splitting given 


59 


17 




31] 


MARKI 


1977 


no splitting given 


21 


18 


l3j 


MARKII 


1979 


splitting given 


24 


19 


m 


MARKI 


1981 


splitting given 


78 



Table 7: Complete data set of measurements of the total hadronic i?-ratio used in this work. Each 
data set is assigned a number to simplify referencing in discussions and figures of this work. Also 
given is information on the name of the experimental collaboration, the year of the publication, 
on whether the splitting of the systematical error in correlated and uncorrelated contributions is 
provided, and on the number of data points. 



4.8 GeV and 6.964 GeV < E < 10.538 GeV. Boughezal et al. |], Kuhn et al. ||, and 



Narison [63] use only one data set from BES [ IT | . Kuhn et al. [35, 36] used the data sets 
of Refs. |l7], |3 from BES covering the energy region 2.6 GeV < E < 4.8 GeV. 

We consider three different selections of data sets to study the dependence of the 
experimental moments on this choice: 

1. The minimal selection contains all data sets necessary to cover the whole energy re- 
gion between 2 and 10.538 GeV without any gaps and keeping only the most accurate 



lg |21|, CrystalBall |2J, CLEO |27j, g6j 



ones. These 8 data sets are from BES |16|, |lj 
and MD1 Q corresponding to the data sets 1, 2, 5, 6, 9, 12, 13, and 14 (see Tab. @ 
for references). 



2. The default selection contains all data sets except for the three ones with the largest 
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uncertainties. It contains 16 data sets and fully includes the minimal selection. It 
contains all data sets except for Mark I and II data sets 16, 17 and 19 from Refs. [33, 
II 

3. The maximal selection contains all 19 data sets. 

We use the default selection as our standard choice for the charm mass analysis, but we 
will also quote results for the other data selections. 

Perturbative QCD region 

Above 10.538 GeV there are no experimental measurements of the total hadronic i?-ratio 
that might be useful for the experimental moments. In this energy region we will therefore 
use perturbative QCD to provide estimates for the charm production ii-ratio. As a penalty 
for not using experimental data we assign a 10% total relative uncertainty to the contri- 
bution of the experimental moments coming from this region, which we then treat like an 
uncorrelated experimental uncertainty for the combination with the moment contribution 



from lower energies. As we see in Sec. |3,3| the energy region above 10.538 GeV contributes 
only (6,0.4,0.03,0.002)% to M^ p 3 4 . In the first moment M^ xp the total contribution is 
about three times larger than the combined statistical and systematical (true) experimental 
uncertainties from the other energy regions. So the 10% penalty we assign to this approach 
represents a subleading component of the final quoted uncertainty. 11 In the second and 
higher moments the contributions from above 10.538 GeV and the corresponding 

uncertainty are negligible compared to the uncertainties from the lower energy regions. 

As the theoretical formula to determine the moment contribution from the perturbative 
QCD region we use the 0(a 3 s ) non-singlet massless quark cross section including charm 
mass corrections up to 0(rn^./s 2 ): 

Rf c (s) = N c Q 2 c R ns (s,m 2 c (^),n f = 4,a: f= \V^)) , (3.1) 

where 

R ns (s,m 2 c (fi),n f = 4,o™ /=4 (/i),^) (3.2) 
— + — 1.52453 - 2.08333L, + I — 

f \l/ \ TT 



! + — + (—) [1.52453 - 2.08333Z S ] + (—) [-11.52034 + 34028^ - 9.56052LJ 

1 \l/ V TT / 

m 2 (u) ( a s /a,\ 2 r , /a s \ 3 r 9 

■|12 — + ^— J [109.167 - 49L S ] + ^— J [634.957 - 799.36LL S + 151.083L, 

+ I^A> f_ 6 + ^£ r_ 22 + 24LJ + (— V [140.855 + 221. 5L S - 5.16667L m - 73L 2 } 

S Z I TT \ TT ) L J 



+ 

S 



+ J [3776.94 - 509.226L S - 174.831L m - 1118.29^ 
+10.3333L^ + 42.1944L m L s + 198.722L^] } , 



11 The situation is quite different if perturbative QCD is used already from 4.8 GeV. Here the contribution 
to the moments M| X 2 P 3 4 from energies above 4.8 GeV amounts to (30, 10, 3, 1)%. For the first moment the 
10% penalty would then represent the largest source of uncertainty and correspond to an uncertainty in 
the charm mass of around 18 MeV. 
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with 



L, = In 



L m = In 



m 2 c {n) 



(3.3) 



•.(J,*/ \ s 

For the computation of the contribution to the experimental moments we determine m c {\fs) 
and a s (y/s) appearing in Eq. ( |3.1| ) using m c (m c ) = 1.3 GeV and a s (mz) = 0.118 as initial 
conditions. 

It is instructive to examine for the moment contributions from ^/s > 10.538 GeV 
terms related to charm production that we do not account for in Eq. (|3.2[ ). In Tab. || the 
relative size with respect to the full first four moments (in percent) of the most important 
neglected contributions are given. In the second column the size of the mass corrections up 
to order m 4 , which we have included in are shown as a reference. The third column 
shows the contributions coming from secondary cc radiation through gluon splitting. The 
fourth column depicts the contributions from the O(a^) singlet corrections (including the 
mass corrections up to order fn^), which one can take as an rough estimate for the actual 
contributions from the charm cut. Finally in the last column we show the size of the Z-boson 
exchange terms integrated from threshold to 10.538 GeV. This contribution represents the 
Z-exchange contribution that is contained in the data, but - by definition - not accounted 
for in the theory moments. We see that at least for the first two moments, the contributions 
neglected are much smaller than the charm mass corrections we have accounted for in the 
nonsinglet production rate, which are already constituting a very small effect. Overall the 
numerical effect on the charm mass of all these contributions is tiny considering the scaling 
Aft'"*. Since we assign a 10% error on the moments' contribution from the energy 
region y/s > 10.538 GeV where we use theory input, our approach to neglect subleading 
effects is justified. 



n 


Mass corrections 


Secondary Radiation 


Singlet 


Z-boson 


1 


0.2 


0.038 


3 x 10" 4 


0.006 


2 


0.01 


9 x 10" 4 


2 x 10~ 5 


0.004 


3 


7 x 10~ 4 


4 x 10~ 5 


2 x 10~ 6 


0.003 


4 


5 x 10" 5 


3 x 10~ 6 


1 x 10~ 7 


0.003 



Table 8: Relative size (in percent) of some subleading contributions for the first four moments 
originating from energies s/s > 10.538 GeV. The second column shows the charm mass corrections 
contained in R^ as a reference. The third column shows the contributions from 0(a 2 s ) secondary 
cc radiation through gluon splitting. The effects from secondary charm radiation accounting all 
energies from threshold to infinity are (0.042, 0.002, 3 x 10 -4 , 7 x 10 _5 )% for the first four moments. 
The fourth column depicts the contribution from the O{o 3 s ) singlet gluon corrections accounting 
for the mass corrections up to order m 4 . Integrating the known singlet corrections from threshold 
to infinity the contributions amount to (0.005, 0.003, 0.002, 0.001)% for the first four moments. The 
last column shows the relative corrections from the Z-boson exchange integrated from threshold to 
10.538 GeV. 



Non-charm background 

Experimentally only the total hadronic cross section is available. Although charm-tagged 
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rate measurements are in principle possible 27] they have not been provided in publications 



On the other hand, they would also exhibit sizable additional uncertainties related to the 
dependence on simulations of the decay of charmed mesons into light quark final states. So 
to obtain the charm production cross section from the data we have to subtract the non- 
charm background using a model based on perturbative QCD related to the production of 
u, d and s quarks. A subtle point is related to the secondary radiation of cc pairs off the u, 
d and s quarks from gluon splitting and to which extent one has to account theoretically 
for the interplay between real and virtual secondary cc radiation which involves infrared 
sensitive terms |34| . Since in this work we define the moments from primary cc production 
(see Eq. (|1.1|)), secondary cc production is formally counted as non-charm background. 
Thus for the model for the non-charm background for ^/s > 2 m c we employ the expression 

Ruds(s) = N c (Ql + 2Q 2 d )R ns (s,0,nf = 3,a¥ = \s)) (3.4) 

/ rif=4 . \ 2 



The second term on the RHS describes the contributions from real and virtual secondary 
cc radiation. The analytic expressions for p R and p v can be found in Eqs. (2) and (6) 
of Ref. ||S5f . We have checked that the numerical impact of real (p R ) and virtual (p v ) 
secondary radiation individually as well as the complete second term on the RHS of 
Eq. ( |3.4p on the moments is negligible, see Tab. ^. We use Eq. fl3.4| ) and fit the non- 
charm background including also data in the region 2 GeV < E < 3.73 GeV via the ansatz 
Rnon-cc(s) = n-ns Ruds (s) > where the constant n ns represents an additional fit parameter. 



3.2 Data Combination 

Combining different experimental measurements of the hadronic cross section one has to 
face several issues: (a) the measurements are given at individual separated energy points, 
(b) the set of measurements from different publications are not equally spaced, cover dif- 
ferent, partly overlapping energy regions and have different statistical and systematical 
uncertainties, (c) the correlations of systematical errors are only known (or provided) for 
the data sets within each publication, (d) there are a number of very precise measurements 
at widely separated energies. 

In this section we discuss the combination of the experimental data from the threshold 
and the data continuum regions between 2 and 10.538 GeV using a method based on a 
fitting procedure used before for determining the hadronic vacuum polarization effects for 
g — 2 [fSlfl . In this work we extend this approach and also account for the subtraction of 
the non-charm background. 

Combination method 

The method uses the combination of data in energy bins (clusters) assuming that the 
R- value within each cluster changes only very little and can thus be well approximated by a 
constant. Thus clusters for energies where R varies rapidly need to be small (in this case the 
experimental measurements are also denser). The R- value in each cluster is then obtained 
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by a x 2 fitting procedure. Since each experimental data set from any publication covers an 
energy range overlapping with at least one other data set, the clusters are chosen such that 
clusters in overlapping regions contain measurements from different data sets. Through 
the fitting procedure correlations are then being communicated among different data sets 
and very accurate individual measurements can inherit their precision into neighbouring 
clusters. Both issues are desirable since the hadronic i?-ratio is a smooth function with 
respect to the sequence of clusters. 

To describe the method we have to set up some notation: 

• All measurements R(E) are distinguished according to the energy E at which they 
have been carried out. 

• Each such energy point having a measurement is written as E^ ,m , where k = 1, . . . , N cx 
refers to the -/V exp data sets, m = 1, . . . , N c i ustev runs over the N c \ uster clusters and 
i = 1, . . . , N k,m assigns the i-th. of the N k,m measurements. 

• Each individual measurement of the -R-ratio is then written as 

R\ E i ) = R i ± <V ± A i > ( 3 - 5 ) 

where R^ ,m is the central value, a^ ,m the combined statistical and uncorrelated sys- 
tematical uncertainty and A k,m the correlated systematical experimental uncertainty. 

• For convenience we define Af^ m = A^ ,m /i?^' m to be the relative systematical corre- 
lated uncertainty. 

As our standard choice concerning the clusters we use 5 different regions each having 
equidistant cluster sizes AE. The regions are as follows: 

(0) non-charm region: has 1 cluster for 2 GeV < E < 3.73 GeV (AE = 1.73 GeV). 

(1) low charm region: has 2 clusters for 3.73 GeV < E < 3.75 GeV (AE = 10 MeV). 

(2) il;(3S) region/threshold region 1: has 20 cluster for 3.75 GeV < E < 3.79 GeV 
(AE = 2 MeV). 

(3) resonance region 2: has 20 cluster for 3.79 GeV < E < 4.55 GeV (AE = 38 MeV). 

(4) continuum region: has 10 cluster for 4.55 GeV < E < 10.538 GeV (AE = 
598.8 MeV). 

We assign to this choice of 52 + 1 clusters the notation (2,20,20,10) and later also examine 
alternative cluster choices demonstrating that the outcome for the moments does within 
errors not depend on them. The cluster in the non-charm region is used to fit for the 
normalization constant n ns of the non-charm background contribution, see Eq. ( |3.4| ). 
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Our standard procedure to determine the central energy E m associated to each cluster 
is just the weighted average of the energies of all measurements falling into cluster m, 



i k,m\<y . / \ k'in\o 

E m = ^ } } • (3-6) 



2 



The corresponding i?-value for the charm cross section that we determine through the fit 
procedure described below is called 

Rm = Rcc{E m ) • (3.7) 

We note that using instead the unweighted average or simply the center of the cluster has 
a negligible effect on the outcome for the moments since the clusters we are employing are 
sufficiently small. 

Fit procedure and x 2 -function 

We determine the charm cross section R C c from a x 2 -function that accounts for the positive 
correlation among the systematical uncertainties A^' m within each experiment k and, at 
the same time, also for the non-charm background. To implement the correlation we 
introduce the auxiliary parameters dk (k = 1, . . . , N exp ) that parametrize the correlated 
deviation from the experimental central values R^' m in units of the correlated systematical 
uncertainty A 4 fc ' m , see Eq. (|3.5|) . In this way we carry out fits to R^' m + dj. A^' m and treat 
the dk as additional auxiliary fit parameters that are constraint to be of order one (one 
standard deviation) by adding the term 

iVexp 

Xcorr({4}) = J2 d l> ( 3 - 8 ) 
fc=l 

to the x 2 -function. To implement the non-charm background we assume that the relative 
energy dependence of the non-charm cross section related to primary production of u, d 
and s quarks is described properly by R u ds given in Eq. ([3.4|). We then parametrize the 
non-charm background cross section by the relation 

Rnon-cc{E) = n ns R u ds(E) , (3.9) 

where the fit parameter n ns is determined mainly from the experimental data in the first 
clusters below 3.73 GeV by adding the term 

XL(-»W) = E E ' ~' ^ * ) "" JWE - » ) , (3.10) 
k=i i=l V a i J 

to the x 2 -function. The complete x 2 -function then has the form 

X 2 ({Rm},n ns ,{d k }) = Xcorr({4}) + {dk}) + X 2 ({i? m },n ns ,{4}), (3.11) 
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where 12 



xl({Rm},n na ,{d k })= (3.12) 

/V cxp N clustc[s N k ' m 

E E E 



■*V -(1 + A/^ d k ) (Rm + n ns RudsiEt' )) 



fc=l m=2 i=l 



k,m 



Note that in our approach the non-charm normalization constant n ns is obtained from a 
combined fit together with the cluster values R m . 

This form of the x 2 " mnc tion is an extended and adapted version of the ones used 
in Refs. [41, 42 1. A special characteristic of the x 2 -functions in Eqs. ( |3.10| ) and ( 3.12] ) 



is that the relative correlated experimental uncertainties Af i ' enter the x -function by 
multiplying the fit value R m rather than the experimental values R k ' m . This leads to a 
non-bilinear dependence of the x 2 -function on the d m and the R m fit parameters and avoids 
spurious solutions where the best fit values for the i? m are located systematically below the 
measurements. Such spurious solutions can arise for data points with substantial positive 
correlation when x 2 -functions with strictly bilinear dependences are employed [41, 42|. 13 



We also note that the implementation of the non-charm background subtraction given 
in Eq. ( p. 12 ) leads to a partial cancellation of systematical uncertainties for the R m best 



fit values for the charm cross section. Moreover, it is interesting to mention that in the 
limit where each cluster contains exactly one measurement (except below threshold, in 
which we always keep one cluster) the x 2 -function decouples, after performing the change 
of variables R' m = R m + n ns R U( i s (E k,m ), into the sum of two independent x 2 -functions, 
one containing data below threshold and depending only on n s and d/~, and another one 
containing data above threshold R k > m (we drop the label i because having only one data 
per cluster it can take only the value 1) and depending only on R' m . After minimizing 
the first x 2 -function one can obtain the best fit values for n s and d k , denoted by ni°^ and 
d^ , respectively. The second \ 2 nas a minimal value of and the best fit parameters read 
B$ = R k ' m /(1 + 4 0) A fe ' m ) - nf ] R uds {E k ^). 

Close to the minimum the x 2 -f unc ti° n of Eq. ( [3.11 ) can be written in the Gaussian 
approximation 



X 2 (fe}) = xlin + £(» -vf ] )Vr\ P] -vf) + 0({p-p^f) , (3.13) 



hi 

where pi = ({R m }, n ns , {c4}) and the superscript (0) indicates the respective best fit value. 
After determination of the correlation matrix Vij by numerically inverting V^" 1 we can 
drop the dependence on the auxiliary variables n nc and d k and obtain the correlation 
matrix of the R m from the i? m -submatrix of Vij which we call , . In order to separate 
uncorrelated statistical and systematical uncertainties from correlated systematical ones 
we compute the complete Vjf, , accounting for all uncertainties and a simpler version of 



We have checked that the effect of using R u ds{E m ) instead of R u d s {E i ' m ) is totally negligible. 
13 We prove in Appendix B the equivalence of a bilinear x 2 -function with fit auxiliary parameters to a 
bilinear x 2 -function without auxiliary fit parameters, but containing the standard correlation matrix. 
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the correlation matrix, V3^, accounting only for uncorrelated uncertainties. The latter is 



obtained from dropping all correlated errors A^" 1 from the x 2 - mnc tion ( [3.11 ) 



11 



The outcome of the fit for the sum of the charm and the non-charm cross section in 
the threshold and the data continuum region using the standard data set explained above 
is shown in Figs. |9|(a)-(f) together with the input data sets. The red line segments connect 
the best fit values and the brown band represents the total uncertainty. The clusters are 
indicated by dashed vertical lines. For completeness we have also given all numerical results 
for the R m values in Appendix A. There we also give results for the minimal and maximal 
data set selections. 

3.3 Experimental Moments 
Narrow resonances 

For the J/ip and V' charmonium contributions to the experimental moments we use the 
narrow width approximation, 

M res = 97rFee (3 14) 

with the input numbers given in Tab. ||. We neglect the tiny uncertainties in the charmo- 
nium masses as their effects are negligible. 

Threshold and data continuum region 

For the determination of the moment contributions from the threshold and the continuum 
region between 3.73 and 10.538 GeV we use the results for the clustered cc cross section 



values R m determined in Sec. 3^ and the trapezoidal rule. Using the relations ds = dE 2 = 
2EdE we thus obtain 

-^clusters t-> r~f 

M , hr+ cont = £ -^(E l+1 -E^ 1 ) + -^L I (E 1 -E ) (3.15) 
i=l E i E o 



lusters + 1 / ITi T-i 
+ E 2n+l l^clustcrs + l ~ &N C 



^clusters ~\~^- 

where Rq and Eq are the R and energy values at the lower boundary of the smallest energy 
cluster, and Rn c1 +i and Em c1 +i are the corresponding values of the upper boundary of the 
highest energy cluster. The values for Rq and i?7v ci +i are obtained from linear extrapolation 
using the R m values of the two closest lying clusters 15 m! and m! + 1 with the formula 

R{E) = ^ m,+1 ~^ m ' (E - E m ,) + R m , . (3.16) 

For the computation of the moment contributions from subintervals within the range be- 
tween 3.73 and 10.538 GeV we also use Eq. ( 3.15|) using corresponding adaptations for the 
boundary values at m = and m = N c i + 1. 

14 In a very good approximation, it can be also obtained by dropping in the rows and columns 

corresponding to dt and inverting the resulting matrix. After that one also drops the raw and column 
corresponding to n nc . We adopt this simplified procedure for our numerical fits. 

15 For Rq we have m' = 1 and for Rn ,+i we have m' = N c i — 1. 
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Whole fit range 



Non-charm region 




(a) 



Low charm region 




3.730 3.735 3.740 3.745 

E (GeV) 

(c) 

Resonance region 




3.750 



3.8 3.9 4.0 




(e) 



J i— 1- 



* i i i » i 



2.0 



2.5 3.0 

E (GeV) 



3.5 



(b) 



(A(3S)/threshold region 




3.75 3.76 3.77 

E (GeV) 

(d) 

Continuum region 



3.78 



3.79 




E (GeV) 



(f) 



Figure 9: Result of the fit for the default selection of data sets. On the top, (a) and (b) show the 
entire fit region and the non-charm region, respectively. In the middle row, (c) illustrates the low 
charm region and (d) the threshold region 1. In the bottom line (e) and (f) depict threshold region 
2 and the data continuum region, respectively. 



Perturbative QCD region 

For the region above 10.538 GeV where we use the perturbative QCD input described in 
Eqs. ( |3.1| ) and ( |3.2| ) for the charm .R-ratio the contribution to the experimental moment is 
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obtained from the defining equation fll.2| ) with a lower integration limit of 10.538 GeV: 



M QCD 



7 x 



ds 



q n+l 



(3.17) 



'(10.538 GeV) 2 s ' 

The variable 7 is an auxiliary variable used to parametrize the 10% uncertainty we assign 
to the perturbative QCD contribution, 



7 



1.0 ± 0.1 



(3.18) 



Correlations 

The experimental moments are obtained from the sum of the resonance, threshold plus 
continuum and perturbative QCD contributions described just above, 



M exp 

n 



M res 
n 



+ M, 



thr+cont 



+ M„ QCD . 



(3.19) 



To determine the uncertainties we account for the errors in the e + e~ widths of J/ip and 
ip' and in the cluster values R m , and for the 10% assigned uncertainty in M„ CD . For the 
evaluation we use the usual error propagation based on a m x fh correlation matrix with 
m = N c i + 3. The correlation matrix of the experimental moments thus has the form 

Ncl+3 /awr p 



a 



exp 



V J V 8R 



( 



R 



(3.20) 



where we have Ri = ({R m },T e + e -(J/ip),r e + e -(ip'),j). The entries of V R in the R m 
subspace are just the entries of the correlation matrix V R obtained from the cluster fit. 
The diagonal entries in the T e + e - subspace are the combined statistical and systematical 
uncertainties of the e + e~ widths and the J7 = 0.1 for Af^ CD , respectively. We treat the 
latter uncertainty as uncorrelated with all other uncertainties. So all non-diagonal entries 
of V R for i or j = N c i + 3 are zero. For the uncertainty of the e + e~ widths we adopt 
a model where the (quadratic) half of the error is uncorrelated and the other (quadratic) 
half is positively correlated among the two narrow resonances and to the R n cluster values 
as well. Thus for the corresponding non-diagonal entries of Vfj with i S {l,N c i} and 
j = {N d + l,N d + 2} we have the entries V$ = 5Tl'+ e _AR m /\/2, and for % = N d + 1 and 
j = N d + 2 we have V R = 5T^ +e _ST^ +e _ /2, where ST^ +e _ are the respective e + e~ width 
total uncertainties for J/ijj and tp', respectively, and AR m the correlated uncertainty of the 
respective cluster R m value. 

The results for the moments showing separately the contributions from the resonances, 
various energy subintervals and their total sum using the defaults data set collection (see 



Sec. 3.1 ) are given in Tab. Using Eq. ( 3.20 ) it is straightforward to compute the corre- 
lation matrix of the moments, and we obtain 

( 0.250 0.167 0.147 0.142 \ 
0.167 0.120 0.107 0.103 
0.147 0.107 0.095 0.092 
0.142 0.103 0.092 0.090 



£<exp 



(3.21) 



7 
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for the total correlation matrix accounting for all correlated and uncorrelated uncertainties. 
We remind the reader that all numbers related to the moment M^ xp are given in units of 
]_q— (n+i) G e v _2n . To quote correlated and uncorrelated uncertainties separately it is also 
useful to show the correlation matrix that is obtained when only uncorrelated uncertainties 
are accounted for. The result it 



y'-fexp 
uc 



/ 0.041 0.035 0.034 0.034 \ 
0.035 0.034 0.034 0.035 
0.034 0.034 0.035 0.036 
0.034 0.035 0.036 0.037 



(3.22) 



These results can be used to carry out combined simultaneous fits to several of the moments. 
This is, however, not attempted in this work. 



n 


Resonances 


3.73-4.8 


4.8-7.25 


7.25 - 10.538 


10.538-oo 


Total 


1 


12.01(17|21) 


3.20(4|6) 


3.47(8|13) 


1.43(2|5) 


1.27(0(13) 


21.38(20|46) 


2 


11.76(18|21) 


1.76(2|3) 


1.13(3|4) 


0.204(3|7) 


0.057(0|6) 


14.91(18(29) 


3 


11.69(19|21) 


0.99(1|2) 


0.390(9|12) 


0.0305(5|10) 


0.0034(0|3) 


13.11(19(25) 


4 


11.77(19|21) 


0.565(8|11) 


0.143(3|4) 


0.00474(9|15) 


2.3(0(2) x 10~ 4 


12.49(19|23) 



Table 9: Result for the experimental moments for the standard selection of data sets. The second 
column collects the contribution from the narrow resonances J/i/j and ip' treated in the narrow width 
approximation. Third to fifth columns are shown only as an illustration, but use the outcome of 
the fit to the entire fit region, 3.73 — 10.538 GeV. The sixth column is obtained using perturbation 
theory, and the last column shows the number for the complete moment. For the moment M^ xp all 
numbers are given in units of 10 - (™ +1 ) GeV~ 2 '\ 



11 


Resonances 


3.73-4.8 


4.8-7.25 


7.25 - 10.538 


10.538-oo 


Total 


1 


12.01(17|21) 


3.11(6|8) 


3.30(9(16) 


1.40(2(6) 


1.27(0|13) 


21.09(22|51) 


2 


11.76(18|21) 


1.73(3|4) 


1.06(3|5) 


0.199(4|9) 


0.057(0|6) 


14.80(19|31) 


3 


11.69(19(21) 


0.98(2|3) 


0.36(1|2) 


0.030(1|1) 


0.0034(0|3) 


13.07(19|26) 


4 


11.77(19|21) 


0.57(1(1) 


0.131(5|5) 


0.0046(1|2) 


2.3(0|2) x 10~ 4 


12.47(19|23) 



Table 10: Results for the experimental moments for the minimal selection of data sets. Conventions 
are as in Tab. ||. 



11 


Resonances 


3.73-4.8 


4.8-7.25 


7.25 - 10.538 


10.538-co 


Total 


1 


12.01(17|21) 


3.16(3|5) 


3.66(6|7) 


1.56(2(4) 


1.27(0|13) 


21.65(19|39) 


2 


11.76(18|21) 


1.75(2|3) 


1.16(2|2) 


0.222(3|5) 


0.057(0|6) 


14.95(18|27) 


3 


11.69(19(21) 


0.98(1|2) 


0.40(1|1) 


0.033(1|1) 


0.0034(0|3) 


13.11(19|24) 


4 


11.77(19|21) 


0.56(1(1) 


0.142(3|2) 


0.0051(1|1) 


2.3(0|2) x 10~ 4 


12.48(19(23) 



Table 11: Results for the experimental moments for the maximal selection of data sets. Conven- 
tions are as in Tab. 0. 
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n 


This work 


Kuhn et al.'07 §6] 


Kuhn et al.'Ol |62] 


Hoang k Jamin'04 [34] 


1 


21.38(20(46) 


21.66(31) 


20.65(84) 


20.77(47(90) 


2 


14.91(18|29) 


14.97(27) 


14.12(80) 


14.05(40(65) 


3 


13.10(19)25) 


13.12(27) 


12.34(79) 


12.20(41(57) 


4 


12.49(19(23) 


12.49(27) 


11.75(79) 


11.58(43(53) 



Table 12: Comparison of our results for experimental moments to those from previous publications 
with a dedicated data analysis. The second column contains our results using the default setting. 
The third and fourth columns were determined from the same collaboration and use data from 
Refs. Jn], [IJJ and jL7|, respectively. The results in the fourth column were also used in the charm 
mass analysis of Ref. 0. The numbers in the fifth column used data from Refs. [[L7J The 
moments in the last two columns were obtained using (less precise) experimental data on the narrow 



resonances |37l that has been updated later 60 



Hoang & Jamin'04 



Kuhn et al.'Ol 



Kuhn et al.'07 



This Work 



20.0 20.5 21.0 21.5 
10 2 xMr P (GeV~ 2 ) 



22.0 



(a) 



Hoang & Jamin'04 



Kuhn et al.'Ol 



13.5 



Kuhn et al.'07 



This Work 



14.0 14.5 
10 3 xM 2 exp (Ge\r 4 ) 

(b) 



15.0 



Hoang & Jamin'04 



Kuhn et al.'Ol 



Kuhn et al.'07 



This Work 



11.5 12.0 12.5 13.0 

10 4 xAf 3 exp (GeV~ 6 ) 

(c) 



13.5 



Hoang & Jamin'04 
' Kuhn et al.'Ol 



Kuhn et al.'07 



This Work 



11.0 11.5 12.0 

10 5 xM 4 exp (GeV~ 8 ) 

(d) 



12.5 



Figure 10: Comparison of determinations of the experimental moments. Blue lines refer to analyses 
in which outdated values for the parameters of the narrow resonances |37|] have been used. 



3.4 Examination 

We conclude this section with an examination of some of the choices and assumptions we 
have implemented for the treatment of the experimental uncertainties. Our defaults choices 
include 

(i) treating one (quadratic) half of J/ip and ip' e + e~ partial width uncertainties as un- 
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correlated and the other half as positively correlated among themselves and with the 
correlated uncertainties of the R m cluster values; 

(ii) treating the entire systematical uncertainties of the i?-ratio measurements as corre- 
lated for the data sets where only total systematical uncertainties were quoted; 

(iii) defining the cluster energies E m through the weighed average of measurement energies 
falling into the clusters, see Eq. fl3,6|); 



(iv) using N c i = 52 + 1 clusters distributed in groups of (2,20,20,10) clusters in the energy 
ranges bounded by (3.73,3.75,3.79,4.55, 10.538) GeV (see Sec. 153) and 



(v) using the default data set collection consisting of all data sets discussed in Sec. 3.1 
except for sets 16, 17 and 19 as defined in Tab. [5]. 



In Tab. 13 alternative correlation models are being studied. The second column shows, 
as a reference, the first four moments with our default settings. The third and fourth 
columns show the moments treating the uncertainties of the J/tp and tp' partial widths being 
minimally correlated with the R m values 16 and completely uncorrelated, respectively. In 
the fifth column we display the moments treating, for data sets 1, 6, 13, 15, 16 and 17, one 
(quadratic) half of systematical uncertainties for the -R-values uncorrelated and the other 
half correlated. In the sixth column all R m values of all data sets are treated as completely 
uncorrelated. We see that the central values depend only weakly on the correlation model 
for those data were the corresponding information is unknown. In particular, for the 
determination of the uncertainties the ignorance about the composition of the systematical 
uncertainties in the i?-values from data sets 1, 6, 13, 15, 16 and 17 is not essential. However, 
for quoting the final uncertainties it is important to account for all (known) correlations 
since treating all errors uncorrelated underestimates the uncertainties considerably. 



In Tab. |lj we examine the impact of modifying the definition of the cluster energy E m 
and of changes to the default cluster distribution (2,20,20,10). In the second column we 
display the resulting moments of the default setting. In the third and fourth columns we 
have shown the moments using for E m simply the mean of the energies and the center of 
the cluster, respectively. The resulting differences to the default definition is an order of 
magnitude smaller than the uncertainties and thus negligible. The fifth, sixth and seventh 
columns display the moments using some alternative cluster distributions. The deviations 
for the default choice illustrated in the table are much smaller than the uncertainties and 
typical for all modifications that satisfy the guidelines for viable cluster definitions we have 



formulated in Sec. |3.2| . This demonstrates that the choice of the cluster distribution does 
not result in a bias for the resulting experimental moments. 

Finally, we also examine the dependence of the moments on the data set collections 



as described in Sec. 3A. In Tabs. [H] and 11 the results for the moments are displayed 



16 We use a modified version of the minimal correlation model. The non-diagonal entries of the correlation 
matrix are filled in with r, R m Min 2 {ATi/Ti, AR m /R m }. Here AI\ and AR m represent the systematical 
uncertainties of the width of the narrow resonance and the R value of the m-th cluster, respectively, and 
i = 1, 2 refer to J/tp and if)' . 
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n 


Default 


Minimal Overlap 


No Correlation 


50% Correlation 


Uncorrelated 


1 


21.38(20|46) 


21.38(20|37) 


21.38(20|30) 


21.24(22|47) 


21.09(28|25) 


2 


14.91(18|29) 


14.91(18(25) 


14.91(18(22) 


14.87(19|30) 


14.84(20|21) 


3 


13.10(19|25) 


13.10(19(22) 


13.10(19(22) 


13.10(19|25) 


13.10(19|21) 


4 


12.49(19|23) 


12.48(19|21) 


12.49(19(21) 


12.49(19|23) 


12.49(19|21) 



Table 13: Dependence on the correlation model. In the second column we show again the results 
for our default set up. In third and fourth column we depict the moments when treating the widths 
of the narrow resonances as minimally correlated and uncorrelated to the values of R m , respectively. 
In the fifth column we treat the quadratic half of the systematical errors of data sets 1, 6, 13, 15, 
16 and 17 as correlated. In the sixth column we show the results when all systematical errors of all 
data sets are treated as uncorrelated. (Note that in the latter case there is still correlation coming 
from the narrow resonances and the contributions from the perturbative QCD region). 



n 


Default 


Regular average 


Middle point 


(2,20,40,10) 


(2,10,20,10) 


(2,20,20,20) 


1 


21.38(20|46) 


21.37(20|46) 


21.40(20|45) 


21.38(20|46) 


21.39(20|46) 


21.34(20|46) 


2 


14.91(18|29) 


14.90(18|29) 


14.91(18|29) 


14.91(18|30) 


14.91(18|29) 


14.88(18|29) 


3 


13.10(19|25) 


13.10(19|25) 


13.11(19|25) 


13.11(19|25) 


13.11(19|25) 


13.08(19|25) 


4 


12.49(19|23) 


12.48(19|23) 


12.49(19|23) 


12.49(19(23) 


12.49(19|23) 


12.47(19|23) 



Table 14: Stability with clustering. In the second column we show the results for our default set 
up. In third and fourth column we depict the resulting moments when for the cluster energy we 
pick the regular average of the energies and the center of the cluster, respectively. In the last three 
columns we change the default clustering to (2, 20, 40, 10), (2, 10, 20, 10) and (2, 20, 20, 20). 



using the minimal and the maximal collections (with default choices for all other settings). 
We see that the differences in the central values to the default collection are about half 
the size of the systematical correlated uncertainties for the first moment M^ xp . For the 
higher moments the differences are much smaller than the uncertainties. Using, instead 
of the default, the minimal and maximal collections affects the uncertainties of M^ xp by 
only about 10%. For the higher moments the differences decrease strongly and basically 
disappear for the fourth moment. Again, the results show that having an increased or 
decreased redundancy in the data set collection only has a minor impact on the final 
numbers for the experimental moments. 

To summarize, we find that modifications to the choices and assumptions that go 
into the combined treatment of the experimental data from different publications and 
experiments lead to changes that are well within the experimental uncertainties we obtain 
from our combination method. We therefore consider these uncertainties as appropriate. 
An instructive comparison of the moments obtained in our analysis to those obtained in 
some previous publications is given in Tab. 12. A graphical illustration of the results is 
shown in Fig. [l(]. 



4. Charm Quark Mass Analysis 

Since it is theoretically most reliable, we use the first moment M\ for our final numerical 
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1.280 - 

1.275 ; 

1.270 

1.265 

1.260 

1.255 



Central value m c (m c ) 



Perturbative Uncertainty 



0.113 0.114 0.115 0.116 0.117 0.118 0.119 

(a) 

Statistical Uncertainty 




0.113 0.114 0.115 0.116 0.117 0.118 0.119 

a s (m z ) 

(b) 

Systematical Uncertainty 




0.00555 



0.113 0.114 0.115 0.116 0.117 0.118 0.119 




0.113 0.114 0.115 0.116 0.117 0.118 0.119 



(c) 



(d) 



Figure 11: Dependence on a s (mz) of the central values of rn c (m c ) (a) and its perturbative (b), 
experimental statistical (c) and experimental systematical (d) errors. The red dots correspond to 
the values we actually calculated. The blue curve shows an interpolation and the red line in (a) 
and (b) corresponds to a linear fit. 



a s (m z ) 


m c (m c ) 


^pcrt 


A st at 


^syst 


m c (3 GeV) 


^pert 


Astat 


^syst 


0.113 


1.255 


0.013 


0.006 


0.013 


1.010 


0.013 


0.006 


0.014 


0.114 


1.259 


0.014 


0.006 


0.013 


1.006 


0.015 


0.006 


0.014 


0.115 


1.263 


0.015 


0.006 


0.013 


1.002 


0.016 


0.006 


0.014 


0.116 


1.267 


0.016 


0.006 


0.013 


0.999 


0.017 


0.006 


0.014 


0.117 


1.271 


0.017 


0.006 


0.013 


0.995 


0.019 


0.006 


0.014 


0.118 


1.275 


0.019 


0.006 


0.013 


0.990 


0.020 


0.006 


0.014 


0.119 


1.280 


0.020 


0.006 


0.013 


0.986 


0.022 


0.006 


0.014 



Table 15: Results for the central values of m c (m c ) and to c (3 GeV) (second and sixth columns) and 
their perturbative (third and seventh columns) , experimental statistical (fourth and eight columns) 
and experimental systematical (fifth and ninth columns) errors. 



charm quark mass analysis. As ingredients for the analysis we use 
(1) the iterative expansion method for the perturbative contribution of the theoretical 
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Central value m ( .(3 GeV) 



Perturbative Uncertainty 




0.113 0.114 0.115 0.116 0.117 0.118 0.119 

(a) 

Statistical Uncertainty 




0.00600 



0.113 0.114 0.115 0.116 0.117 0.118 0.119 

a£rnz) 



(c) 




0.014 



0.113 0.114 0.115 0.116 0.117 0.118 0.119 

a s (m z ) 

(b) 

Systematical Uncertainty 




0.113 0.114 0.115 0.116 0.117 0.118 0.119 



(d) 



Figure 12: Dependence on a s (mz) of the central values of m c (3GeV) (a) and its perturbative 
(b), experimental statistical (c) and experimental systematical (d) errors. The red dots correspond 
to the values we actually calculated. The blue curve shows an interpolation and the red line in (a) 
and (b) corresponds to a linear fit. 



moment at O(a^), see Eq. (^TJ), 
(2) the gluon condensate correction with its Wilson coefficient determined at 0{a s ) as 



described in Sec. 2.2 



(3) the first experimental moment 



M° xp = 0.2138 ± 0.0020 st at ± 0.0046 ; 
using our default settings as discussed in Sec. [3|. 



syst 



(4.1) 



One important source of uncertainty we have not yet discussed is the value of the strong 
MS coupling a s . Since in the recent literature 67, 68, 69. 70 1 a s determinations with 



a spread larger than the current world average [71] have been obtained, we carry out 
our numerical analysis for values of a s (mz) between 0.113 and 0.119. 17 The outcome of 



lr As our default we use ctl 1 (mz) as the input, use the four-loop QCD beta-function for the renormal- 
ization group evolution and three- loop matching conditions to the n/ = 4 theory at pi = 4.2 GeV. 
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our analysis is shown in Tab. 15. In Fig. 11 the central values, perturbative, statistical 



and systematical uncertainties are displayed graphically. For the central value and the 
perturbative uncertainty, which show a significant dependence on a s , we can present a 
linear fit. For the statistical and systematical uncertainties the variation with a s is smaller 
than 1 MeV and we only quote constant values. We thus obtain 

m c (m c ) = (0.788 + 4.13 x a s {m z )) ± (0.006) stat ± (0.013) syst (4.2) 
± (-0.127 + 1.23 x a s (m z )) pcrt ± (0.002) (GG) , 

m c (3 GeV) = (1.462 -4.00 x a s {m z )) ± (0.006) stat ± (0.014) syst (4.3) 
± (-0.148 + 1.43 x a s (m z )) pcvt ± (0.002) (GG) . 



Taking as an input 



a s (m z ) = 0.1184 ± 0.0021 , (4.4) 



which is the current world Bethke average [71] with a tripled uncertainty we obtain 



m c (m c ) = 1.277 ± (0.006) stat ± (0.013) syst ± (0.019) pert ± (0.009) Qa ± (0.002) (GG) , 

m c (3 GeV) = 0.998 ± (0.006) stat ± (0.014) syst ± (0.022) pert ± (0.010) a , ± (0.002) {GG) , 

(4.5) 

which represents, together with Eq. ( [4.2[ ), our final numerical result for the MS charm 
mass. Our result is in good agreement with other recent determinations of m c (m c ). A 



summary of the numerical results is shown in Tab. 16 and in Fig. |13|. Compared to the 
analysis carried out in Refs. |?5], |3(| our experimental uncertainty is larger by 5 MeV 
and our perturbative uncertainty is a factor of 10 larger. Compared to Refs. || [| the 
discrepancy in the perturbative error estimate is even larger. We note that the almost 



head on agreement in the central value of Refs. [35, ^6[ is accidental, since the central 
value we obtain is based on a different perturbative expansion (iterative expansion in our 
work versus fixed-order expansion in theirs), on the way how the central value is defined 
(average of the lowest and highest values in our work versus value for fx = 3 GeV) and on 
a different central value for the experimental moment (M^ xp = 0.2138 in our work versus 
M° xp = 0.2166 in theirs). Whereas our definition of the central value makes m c smaller, 
our smaller value for M° xp makes it bigger by almost the same amount. 



5. Conclusions and Final Thoughts 

In this work we have used state-of-the-art 0{ofi s ) input from perturbative QCD to deter- 
mine the MS charm quark mass from relativistic (n = 1) charmonium sum rules using 
experimental data on the total hadronic cross section in e + e~ annihilation. The main 
aims were (i) to carefully reexamine perturbative uncertainties in the charm mass extrac- 
tions from the moments of the charm vector current correlator and (i) to fully exploit the 
available experimental data on the hadronic cross section. 
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'mc(' m c) ot s (mz) used m c {m c ) aa ^ mz ^ °- 1180 


This work 

Chetyrkin et al. [35] 
Boughezal et al. || 
Hoang & Jamin [34] 
Bodenstein et al. |3l|] 
Narison 63] 
Allison et al. @ 
McNeile et al. [||] 


1.277 ±0.026 0.1184 ±0.0021 1.275 ±0.023 
1.279 ±0.013 0.1189 ±0.0020 1.277 ±0.012 
1.295 ±0.015 0.1182 ±0.0027 
1.29 ±0.07 0.1180 ± 0.0030 1.29 ± 0.07 
1.319 ±0.026 0.1213 ±0.0014 1.295 ±0.026 
1.260 ±0.018 0.1191 ±0.0027 
1.268 ±0.009 0.1174 ±0.0012 
1.273 ±0.006 0.1183 ±0.0007 



Table 16: Comparison of recent determinations of m c (rn c ). In the second column we show the 
final result as quoted in the publications, and in the third the value of ot s (mz) used in the analysis. 
In the fourth column, when possible, we extrapolate the results to a s (mz) = 0.1180. For the result 
of Ref. |n| we perform a linear extrapolation using the two values for a s (mz) = 0.1189 and 0.1213 
quoted in the paper. The last two rows correspond to lattice studies, where the value of a s (mz) 
was simultaneously determined from the analyses. 



We carried out this work having in mind recent 0(a s ) sum rule analyses [||, ^, 35, 



36], where perturbative errors of 2 MeV or smaller were quoted. Moreover in Refs. [35 
] an experimental uncertainty of 9 MeV was quoted. Given these numbers we found 
it appropriate to reexamine their analysis. In their work the perturbative uncertainty 
estimate was achieved using a specific choice to arrange the perturbative expansion and 
by setting the renormalization scales in a s (fi a ) and of the MS charm quark mass (// m ) 
equal. We demonstrated that this results in an accidental cancellation of /i a and \i m 
scale variations that is not observed in other alternative ways to treat the perturbative 
expansion. Moreover, concerning the experimental input their work relied on perturbative 
QCD predictions instead of available data for energies ^fs > 4.8 GeV which results in a 
rather strong dependence of the final experimental uncertainty on the model-uncertainty 
assigned to the theory input. 

Concerning a proper assessment of perturbative and experimental uncertainties we 
established in our work the following improvements: 

1. We demonstrated that for achieving an estimate of perturbative uncertainties based 
on scale variations that is independent of the perturbative expansion method one 
needs to vary [x a and [x m independently, albeit with ranges that avoid large logs. As 
a result the perturbative uncertainty estimates using different ways to carry out the 
expansion in a s become equivalent, which is not the case for [i a = [i m . 



Using a data clustering method similar to Refs. [41, 42, 43] we combined available data 
on the total e + e~ hadronic cross section from many different experiments covering 
energies up to ^fs = 10.538 GeV to fully exploit the existing experimental information 
for the experimental moments. This avoids a strong dependence of the experimental 
moments on ad-hoc assumptions on the "experimental" model uncertainty being 
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Hoang & Jamin'04 
Chetyrkin et al.'09 

Boughezal et al.'06 
This Work 




1.22 


1.24 1.26 1.28 1.30 1.32 1.34 

m c (m c ) 


1.36 



Figure 13: Comparison of recent determinations of m c (m c ). Red corresponds to our result, black 
and gray correspond to Oaf) and 0{a 2 s ) charmonimum sum rules analyses, respectively, green labels 
other kind of sum rule analyses (weighted finite energy sum rules J6l[ | and ratios of Q 2 -dependent 
moments JS3|), and blue stands for lattice simulations. 



associated to the QCD theory input used for the experimental moments. As a result 
we were also able to quantify the correlation between different experimental moments. 

Using a s (m z ) as an unspecified variable and the theoretically more reliable first moment 
Mi for the fits we have obtained 

m c (m c ) = (0.788 + 4.13 x a s (m z )) ± (0.006) stat ± (0.013) syst (5.1) 
±(-0.127 + 1.23 x a s (m z )) pcit ± (0.009) as ± (0.002) (GG> , 

m c (3 GeV) = (1.462 -4.00 x a s {m z )) ± (0.006) stat ± (0.014) syst (5.2) 
± (-0.148 + 1.43 x a s (m z )) pcit ± (0.002) (GG> , 

for the MS charm mass. At the level of uncertainties obtained in our work excellent 
convergence of perturbation theory was observed. Adopting a s (m z ) = 0.1184 ± 0.0021 we 
then obtain 

m c (m c ) = 1.277 ± (0.006) stat ± (0.013) syst ± (0.019) pert ± (0.009) aa ± (0.002) {GG) , 

m c (3 GeV) = 0.998 ± (0.006) stat ± (0.014) syst ± (0.022) pert ± (0.010) a , ± (0.002) (GG) . 

(5.3) 

Our perturbative error of 20 MeV is a factor of ten larger, and the experimental uncertainty 
of 14 MeV is by 5 MeV larger than in the most recent analysis of Ref. [36]. For estimating 
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the perturbative error a range of scale variations between m c (m c ) and 4 GeV was employed. 
Adding all uncertainties quadratically we obtain 

m c (m c ) = 1.277 ± 0.026 GeV, 
m c (3 GeV) = 0.998 ± 0.029 GeV, (5.4) 



giving an uncertainty that is twice the size of the one obtained in Refs. [^g, 36]. This 



difference arises mainly from the more appropriate estimate of perturbative uncertainties 



obtained in this work. The agreement in the central value with Refs. |35|, 36 1 is a coin- 
cidence, arising from an interplay of differences in the experimental input and theoretical 
treatment. 

As a final thought one might ask which further improvements might be possible in the 
future. As can be seen from Tab. |9|, from the experimental side the biggest improvement 
could be made from more accurate measurements of the J/ip and tp' electronic partial 
widths. The current relative uncertainties are 2.5% and 2.4%, respectively. Here some 
improvement might be conceivable with dedicated measurements. On theory side, viewing 
the uncertainties, it is not unreasonable to assume that the computation of 0{a^) moments 
of the vector current correlator might further reduce the perturbative error below the level 
of 20 MeV. Using the MS scheme for the charm mass the OPE states that the remaining 
perturbative infrared renormalon ambiguity is of order AQ CD /mj? ~ 0(5-15 MeV). This 



expectation has also been confirmed by explicit bubble chain (large- (3q) computations [72] 
and indicates that a further reduction of the perturbative uncertainty is not excluded. 

However, to throw in some words of caution, at the level of the present perturbative 
uncertainties one should also remind oneself about possible loopholes still left in the char- 
monium sum rule method. An issue we would like to mention concerns the separation of 
charm and non-charm hadronic production rates needed to carry out the charmonium sum 
rule. On the theory side the issue is conceptually subtle due to the singlet and secondary 
radiation contributions which arise at O(af) and 0{a 2 s ), respectively. In this work (as 
well as in Ref. |]35[] ) both contributions have been considered as non-charm although they 
contain terms belonging to the cc final state. This treatment might be justified since the 
size of the corresponding terms is quite small (see Tab. ||) and since it is the common ap- 
proach to determine the experimental charm production rate in the continuum region by 
subtracting theoretical results (or models) for the non-charm rate from the measured total 
hadronic iZ-ratio. In our method to determine the experimental moments this subtraction 
involves a normalization constant multiplying the theoretical non-charm i?-ratio that is 
fitted within our clustering method as well accounting for data below and above the charm 



threshold. The result (see Eq. ( |A.2|) ) reveals a disparity of 4% between the theoretical 
non-charm i?-ratio and the data. Setting, in contrast, the normalization constant to unity 
results in a shift in the charm mass by —15 MeV. 18 Since this shift is compatible with 



18 The same disparity was found in Ref. | pq |. In their analysis the corresponding effect is —5 MeV since 
experimental data were used in the experimental moments only for energies ^/s < 4.8 MeV. In this approach, 
however, the moments are strongly dependent on the uncertainty one assigns to the theory input used for 
v/i > 4.8 MeV. 
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the overall systematical uncertainty in the experimental data we have not treated it as an 
additional source of uncertainty. On the other hand, the size of the shift could also be 
considered as an inherent conceptual uncertainty related to separating the charm from the 
non-charm .R-ratio, which is based on theory considerations rather than on experimental 
methods and which apparently cannot be improved simply by higher order perturbative 
computations. As an alternative, one might avoid the separation of the charm and the non- 
charm contributions altogether and use the total hadronic cross section for the charm mass 
fits. Apart from the shift mentioned above such an approach would, however, also lead to 
a substantially larger dependence on the uncertainties in a s . Given these considerations 
we believe that a substantial reduction of the uncertainties also relies on a resolution of 
the disparity mentioned above. This might certainly involve more precise measurements in 
the charm threshold and below-threshold regions, but also some deeper conceptual insight. 
Until then a substantial reduction of the uncertainties shown in Eq. ( |5.4| ) appears hard to 
achieve without imposing ad-hoc assumptions. 
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Appendix A: Results of the Fit Procedure 

In this appendix we present in some more detail the numerical results of our fit procedure. 
In Tabs. |A.1| , [A. 2 and [A. 3] , the results for the cluster energies and the cluster charmed 



i?-values are shown for the standard, minimal and maximal selection of data sets, respec- 
tively, using our default setting for the correlations. We use the results for the standard 
data set selection for our final charm mass analysis. The numbers in the parentheses corre- 
spond to the statistical and systematical errors. The correlation matrices for the -R-values 
is available, but cannot be displayed due to lack of space. They can be obtained by the 
authors on request. For the three data selections, the fit gives the following minimal x 2 
per degree of freedom, 

^standard gg X minimal ^ gg --^maximal ^ g-^ 

dof dof dof 
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Figure 14: Comparison of the results for the experimental moments using the three data selections. 



and the following normalization constants for the non-charm background 



, standard 



'IK 



, maximal 



>K 



1.039 ± 0.003 stat ± 0.012 
1.023 ± 0.003 sfc , t ± 0.011 



syst j 
syst- 



, minimal 



I'K 



1.029 ± 0.003 stat ± 0.015 



-•syst ; 



(A.2) 



Since the minimal x 2 /dof values are not close to unity, one has to conclude that fit quality 
is not really very good. This is not at all visible from the agreement of the fit and the data 
for the total cross section (see Figs. and thus might be related to the disparity between 
the fits of charm versus non-charm production rates described in Sec. ||. 

In Eq. (A. 3) we show for the correlation matrices of the first four experimental moments 
for the minimal and the maximal data set selection. The results for our standard selection 
are given in Eqs. ( 3.2l| ) and ( |3.22| ). All numbers are related to moments M„ xp normalized to 
units of 10~( n+1 ) GeV _2n . We show the results accounting for the full set of correlated and 
uncorrelated uncertainties and the correlation matrices accounting only for uncorrelated 
systematical and statistical uncertainties (subscript uc). 
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3.769 


1.31(8|3) 


3.771 


1.39(6|4) 


3.773 


1.42(2 4) 


3.775 


1.22(8|3) 


3.777 


1.29(8|4) 


3.78 


1.13(7|3) 


3.781 


1.06(5|3) 


3.783 


1.00(7|3) 


3.785 


0.63(11|3) 


3.787 


0.64(7|3) 


3.79 


0.41(5|2) 


3.808 


0.11(4|2) 


3.846 


0.10(5|2) 


3.883 


0.24(7|2) 


3.928 


0.63(8|2) 


3.967 


0.94(3|2) 


4.002 


1.30(3|3) 


4.033 


2.18(3|4) 


4.069 


1.74(4|3) 


4.117 


1.72(4|3) 


4.156 


1.68(1|4) 


4.191 


1.54(3|3) 


4.23 


0.91(5|2) 


4.261 


0.71(2|2) 


4.307 


0.89(6|2) 


4.346 


1.18(6|2) 


4.382 


1.63(7)3) 


4.416 


1.76(4)3) 


4.452 


1.53(7|3) 


4.492 


1.42(6|3) 


4.529 


1.29(10|3) 


4.715 


1.51(3|3) 


5.326 


1.39(7|7) 


6.006 


1.33(6|7) 


6.596 


1.35(7|7) 


7.202 


1.39(2|4) 


7.852 


1.52(8|5) 


8.417 


1.38(3|5) 


9.04 


1.43(5|6) 


9.54 


1.35(3|5) 


10.327 


1.37(1|5) 



Table A.l: Best fit values for the standard selection of data sets. The eneri 
measured in GeV, and for R the first number in brackets is the statistical error 
systematical one. 



ry of the cluster is 
and the second the 



£rexp 



/ 0.304 0.194 0.165 0.157 \ 
0.194 0.133 0.115 0.109 
0.165 0.115 0.100 0.096 

\ 0.157 0.109 0.096 0.092 ) 

( 0.193 0.138 0.125 0.122 \ 
0.138 0.108 0.099 0.097 
0.125 0.099 0.091 0.090 
0.122 0.097 0.090 0.088 



/-yexp 



/ 0.047 0.038 0.035 0.035 \ 
0.038 0.035 0.035 0.035 
0.035 0.035 0.035 0.036 
0.035 0.035 0.036 0.037 

( 0.037 0.034 0.033 0.034 \ 
0.034 0.033 0.034 0.035 
0.033 0.034 0.035 0.036 
0.034 0.035 0.036 0.037 



(A.3) 



Appendix B: On the Equivalence of x 2 -Functions 



In this appendix we demonstrate that a x 2 -function in which the auxiliary fit parameters dk, 
which describe the correlated deviation off the experimental central value within experiment 
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E 


R 


E 


R 


E 


R 


E 


R 


3.736 


-0.07(6|2) 


3.749 


0.22(6|2) 


3.751 


0.38(9|3) 


3.753 


0.38(7|3) 


3.755 


0.40(6|3) 


3.757 


0.68(11|3) 


3.759 


0.74(7|3) 


3.761 


0.82(7|4) 


3.763 


1.01(8|4) 


3.765 


1.21(9|4) 


3.767 


1.49(12|4) 


3.769 


1.27(8|4) 


3.772 


1.37(6|4) 


3.773 


1.38(8|4) 


3.775 


1.21(8|4) 


3.777 


1.26(8|4) 


3.78 


1.11(7|3) 


3.781 


1.06(5|4) 


3.783 


0.97(7|4) 


3.785 


0.62(11|3) 


3.787 


0.66(8|3) 


3.79 


0.40(5|3) 


3.808 


0.11(5|2) 


3.845 


0.12(5|2) 


3.879 


0.18(9 2) 


3.935 


0.85(14 3) 


3.969 


1.01(4 3) 


4.002 


1.39(3|4) 


4.032 


2.25(5|5) 


4.066 


1.94(4|5) 


4.118 


1.78(5|5) 


4.157 


1.77(1|5) 


4.191 


1.62(3|5) 


4.232 


0.99(8|3) 


4.261 


0.76(2|3) 


4.309 


0.83(10|3) 


4.35 


1.22(11|3) 


4.385 


1.28(14|4) 


4.415 


1.62(10|4) 


4.45 


1.54(11|4) 


4.492 


1.30(13|4) 


4.529 


1.03(13|3) 


4.716 


1.25(8|3) 


5.378 


1.40(8|8) 


6.008 


1.31(7|8) 


6.622 


1.30(8|8) 


7.206 


1.37(2|6) 


7.856 


1.47(8|7) 


8.417 


1.36(3|6) 


9.037 


1.40(5)7) 


9.544 


1.31(3|6) 


10.252 


1.34(3|6) 



Table A. 2: Best fit values for the minimal selection of data sets. Conventions are as in Tab. lA.lL 



E 


R 


E 


R 


E 


R 


E 


R 


3.736 


0.01(6|2) 


3.749 


0.29(5|2) 


3.751 


0.48(8)3) 


3.753 


0.40(7|2) 


3.755 


0.42(6|3) 


3.757 


0.80(10|3) 


3.759 


0.76(7|2) 


3.761 


0.85(6|3) 


3.763 


1.06(8|3) 


3.765 


1.24(8|3) 


3.767 


1.51(12|4) 


3.769 


1.31(7|3) 


3.771 


1.39(6|4) 


3.773 


1.41(2|4) 


3.775 


1.21(8|3) 


3.777 


1.28(8|4) 


3.78 


1.14(7|3) 


3.781 


1.06(5|3) 


3.783 


1.00(7|3) 


3.785 


0.61(9|3) 


3.787 


0.64(7|3) 


3.79 


0.41(5|2) 


3.808 


0.13(4|2) 


3.846 


0.10(4|2) 


3.884 


0.26(5|2) 


3.927 


0.66(7|2) 


3.967 


0.98(3|2) 


4.002 


1.30(2|2) 


4.033 


2.18(3|3) 


4.07 


1.74(3|3) 


4.117 


1.75(4|3) 


4.156 


1.70(1|4) 


4.191 


1.54(3|3) 


4.23 


0.91(5|2) 


4.261 


0.74(2|2) 


4.307 


0.87(6|2) 


4.346 


1.15(5|2) 


4.382 


1.55(6|3) 


4.416 


1.80(4|3) 


4.452 


1.52(6|3) 


4.492 


1.33(5|2) 


4.53 


1.19(8|2) 


4.722 


1.42(3|2) 


5.36 


1.39(6|4) 


6.018 


1.47(4|3) 


6.608 


1.66(4|3) 


7.202 


1.54(2|3) 


7.851 


1.63(8|5) 


8.417 


1.50(3|4) 


9.04 


1.55(5|5) 


9.54 


1.47(3|4) 


10.327 


1.48(1|4) 



Table A. 3: Best fit values for the maximal selection of data sets. Conventions are as in Tab. A.l 
k, multiplies only the experimental systematical uncertainties, 19 



X 



k 




R,: ■ + d k A • ' 



Rr 



k,m 



(B.l) 



For the discussion in this appendix we do for simplicity of the presentation not account for fits of the 
non-charm contribution. It is, however, straightforward to generalize the presentation to this case. 
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is mathematically equivalent to the well known % 2 -function written solely in terms of the 
fit parameters R m and a correlation matrix, 



i,j,k,m,n 



(B.2) 



where there is no correlation among the different experiments k and k', and the correlation 
matrix within one experiment k has the form 



(Tr\mn k,m2 x rmn , \k,m\k,n 

{Vk)ij =o i dijd + A t Aj . 
The matrix Eq. QB.3]) can be inverted analytically giving 



(Vk )ij 



l\mn _ u iJ 



1 A*'™ A^' n 

2\ / fc.mxo / k,n\ r . 



(°i ) 1 | \ A /fc ( a i ) (^i 

/ \? 



(B.3) 



(B.4) 



where we have defined the mean statistical error and the statistical average of the system- 
atical error within one experiment as follows: 



(1 \ 1 2 (A fe ' m ) 2 

i,m y^i ) J i,m ) 



(B.5) 



In case of a sizable positive correlation between measurements in the same experiment 
(such that A^' m /R^' m is sizable and approximately constant) is it known that the form of 
the x 2 -function in Eq. (|B.1| ) leads to best fit values R m that are systematically below the 
measurements |2l]] . Our proof demonstrates that the standard x 2 -function of Eq. ( BJ2 ) has 
the same property for correlation matrices with the form of Eq. (|B.3|). This motivates to 



use the so-called minimal-overlap correlation model where the second term on the RHS of 
Eq. ( |B.2[ ) is replaced by min 2 (A^' m , A k - )Tl ). In general, within the minimal-overlap model, 
the correlations are sufficiently reduced such that the unphysical effect described above 
does not arise. 

We proceed by showing that one can "integrate out" auxiliary parameters dk in 
Eq. ( |B.1| ) obtaining a new function x 2 (i?*) which yields the same results for the best fit for 
the R m , as long as one works in the Gaussian approximation. The minimum of x 2 {Ri, dk) 
is located at the best fit values (indicated by superscripts (0)) defined by the conditions 



dx 2 



dRi 



dx 2 



r: 



(0) .(0) 



Odi 



0. 



R (0) d {0) 



(B.6) 



To invert the matrix of second derivatives we proceed in blocks 



dRidRj dRiddk 
2 X 2 d 2 x 2 
V ddidRj ddidb k ) s 



b'rnk &kr 



Sir 



(B.7) 



R,b 
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where Cij and Cvky* 3X6! 

-^cluster x -^cluster and iV ex p x iV ex p square matrices, respectively, and 
bj r is a ^duster x -^exp rectangular matrix. We find the following four matrix relations 



E 



2^2 



d 2 x 2 



dRiddk 



dRiddk 



+E 



<9 2 x 2 



Combining Eqs. 



ddiddk 
and 



ddiddh 



= 0, (B. 
<5/ r • (B. 



) we find the inverse of the upper left i?-block, 



where 



(c-% = 2(V-% 

d 2 x 2 



d 2 x 2 



dRjdRi 



E 

kl 



d 2 x 2 



ddidRi 



d 2 x 2 



ddidd k 



-i 



d 2 x 2 



dd k dR, 



d 2 x 2 



stands for the (I, k) element of the inverse matrix of 



(B.10) 
(and not 



ddiddk J ddiddk 
the inverse of the element). Combining Eqs. ( |B.8[ b) and ( |B.9| b) one can obtain relations 

for the a and b submatrices: 

-l 



3 r 



E 



d 2 x 2 



dRjdRi 



d 2 x 2 



6lr = 



d 2 x 2 



ddi ddk 



E 



ddkdR, 

d 2 x 2 

ddidRj 



(Ikr 



i -J 



2, ,2 "1-1 



d 2 x 



dRjdRi 



2, ,2 



d X 



ddidRi 



CLkr > 



d 2 x 2 



E 



5 2 x 2 



dRidRj ~ 31 " ^ ddkdR, 

J k,Ls 



d 2 x 2 



dRmdR.s 



-i 



d 2 x 2 

dd t dR ; 



where again 



d 2 x 2 



stands for the (J, i) element of the inverse matrix of 



(B.ll) 

(B.12) 
(B.13) 

d 2 x 2 

dRjdRi ' 



dRjdRi 

We "integrate out" the auxiliary parameters dk by substituting their minimum conditions 
d k °^ = d~k(Ri°^) (which is analogous to using the equation of motion when integrating out 
heavy particles): 



dx 2 



ddi 



(B.14) 



d k =d k (Ri) 



Their first derivatives with respect to Rj read 



d dx 2 



dRj ddi 



d 2 x 2 



dm(Rk) 



ddidRj 



ddj(R n ) 
dRj 



E 



. +E 

dm(Rk) I 

2 i-l 



d 2 x 2 



d 2 x 2 

ddidd, 



ddiddi 

d 2 x 2 



ddj(R rt 
dm(Rk) dR i 



The minimum of x (Ri) is indeed located at R^ because 



dm (Rk ) 

(0) 



ddidRj 



dm(Rk) 



(B.15) 



dx 2 _ \- d X 2 
dRi V dd k 



dd k dx 2 



dk-dk(Ri) 



dRi dRj 



dx 2 



dk-dk(Ri) 



dRi 



dk-dk(Ri) 



(B.16) 
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and because the first term vanishes by the condition in Eq. ( B.14 ). When evaluating 
Eq. ( B.16 ) for B4 = R^ it vanishes because of Eq. ( B.6 ). Finally, let us calculate the 



inverse correlation matrix: 



d 2 x 2 



dRidRj 



d 2 x 2 



d 2 x 2 



R k dRidRj 



+ £ 



d 2 x 2 



ddk(Ri) 



(o) D (o) dRj 



dm ,R 



Rk 



dRidRj 



E 

k,l 



d 2 x 2 

dd k ddi 



d 2 x 2 



d (o) 0(0) ddidRi 



d 2 x 2 



d (o) ii? (oj ddkdRj 



AO) R (o) 



(B.17) 



which agrees with Eq. (B.10). 

We can now apply the previous results to Eq. 

2" 



X 



£ 



4 + E 



R k,m + dkA k,m_ Rm 



k,m 



(B.18) 



£4 



1 + 



+^E^^ + Ef* m 



k,m 



The equation of motion for reads 
d k (R) = 



£ 



(B.19) 



1 + 



<*>2 



This renders for x 2 the form 



x 2 



E 



E 



4(#)£ 



{Rf m - i? m )A* 



fc,m 2 



+E 



k,m 



Rr, 



k,m 



(7; 



(B.20) 



E - 



k,m 



R r , 



k,m 



0"; 



1 + 



I A 2 



E 

i,j,m,n 



(R^ m -R m )A^ m {Rj' n — Rn)A 



. k,n 



a 



k,m 2 



(7 



k,n2 



which reproduces Eq. (B.2), as we wanted to demonstrate. 

We conclude this appendix presenting an alternative way to write Eq. ( |3.11| ) after 
using the equations of motion for dk- We again concentrate on the simpler case without 
subtraction of the non-charm contribution: 

"i^ m -(i+4A/r)^ 12 



X' 



E 4 + E 



= E 

k 

+ E 



!+E 



itu — it 



i.m \ k 
2 ' 



4 - 2d k 



Rl 



Rm ) A/?' m i? r , 



K ) 



7(1 



(Til 



(B.21) 
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The EOM for du now reads 



1 + ^(AflTR» 



cr. 



^ K K k ~ timj^h K ™ 



K ) 



which upon insertion into Eq. (B.21) renders 

2 



X 



1 i,m 


- n 'j! c Jrt m 


i,m 


1. i,m 







Rm ) A/^' m i? m 



K ) 



E 



cr 



One can rewrite Eq. (B.23) in the matrix form 



(B.22) 



(B.23) 



X 



k I i,j,m,n 



ft" 



.nmn OjjO 



i?,„) \f k ' m R m [R{ n -R n )Ml n B„ 



R 



K ) 



i + E 



(B.24) 



In Eq. flB~2l ) one can interpret the second term of the inverse correlation matrix, as a non- 
linear correlation among the measurements, where the correlation matrix itself depends on 
the value of the fit parameters. The total inverse correlation matrix is block diagonal, and 
the blocks correspond to VC ■ 
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